module weno_1d_mod
  use recon_type_mod
  use recon_math_mod
  use poly_utils_mod
  use smooth_indicators_mod

  implicit none

contains

  recursive subroutine weno_1d_init(this, nd, sw, swx, swy, swz, xc, yc, zc, dx, dy, dz, ic, jc, kc, is, ie, js, je, ks, ke, mask, id)

    class(recon_type), intent(inout) :: this
    integer , intent(in), optional :: nd
    integer , intent(in), optional :: sw
    integer , intent(in), optional :: swx
    integer , intent(in), optional :: swy
    integer , intent(in), optional :: swz
    real(16), intent(in), optional :: xc(:)
    real(16), intent(in), optional :: yc(:)
    real(16), intent(in), optional :: zc(:)
    real(16), intent(in), optional :: dx
    real(16), intent(in), optional :: dy
    real(16), intent(in), optional :: dz
    integer , intent(in), optional :: ic
    integer , intent(in), optional :: jc
    integer , intent(in), optional :: kc
    integer , intent(in), optional :: is
    integer , intent(in), optional :: ie
    integer , intent(in), optional :: js
    integer , intent(in), optional :: je
    integer , intent(in), optional :: ks
    integer , intent(in), optional :: ke
    integer , intent(in), optional :: mask(:,:,:) ! Cell mask
    integer , intent(in), optional :: id

    integer i, j, k, m, n, sub_ie, sub_je, ingb, isub
    integer max_ngb

    call this%clear()
    
    if(nd/=1)then
      stop 'weno_1d_init support only nd==1'
    endif
    
    this%sw      = sw
    this%swx     = sw
    this%swy     = 1
    this%swz     = 1
    this%nd      = nd
    this%nc      = sw**nd
    this%max_ngb = 3**nd-1
    this%rw      = ( sw - 1 ) / 2
    
    max_ngb = this%max_ngb
    
    this%rwL = this%rw
    this%rwR = this%rw
    
    this%dx = dx
    if(present(dy))then
      this%dy = dy
    else
      this%dy = 1
    endif
    if(present(dz))then
      this%dz = dz
    else
      this%dz = 1
    endif
    
    this%eps = 1.e-15
    
    if (present(id)) this%id = id

    if(present(is))then; this%is = is; else; this%is = 1       ; endif
    if(present(ie))then; this%ie = ie; else; this%ie = this%swx; endif
    if(present(js))then; this%js = js; else; this%js = 1       ; endif
    if(present(je))then; this%je = je; else; this%je = 1       ; endif
    if(present(ks))then; this%ks = ks; else; this%ks = 1       ; endif
    if(present(ke))then; this%ke = ke; else; this%ke = 1       ; endif

    allocate(this%cell_mask(this%is:this%ie,1,1))
    allocate(this%poly_mask(sw,1,1))

    if (present(mask)) then
      this%cell_mask = mask
    else
      this%cell_mask = 1
    end if

    ! Set polynomial mask from cell mask.
    if (this%cell_mask(this%is,1,1) == 0) then
      this%poly_mask = this%cell_mask(this%ie:this%is:-1,:,:)
    else
      this%poly_mask = this%cell_mask
    end if

    ! Only check trouble cell on stencil.
    if (this%id == 0) then
      allocate(this%cell_ngb(nd,max_ngb,sw,1,1))
      allocate(this%cell_ngb_size(sw,sw,1))
      k = 1
      j = 1
      do i = this%is, this%ie
        if (this%cell_mask(i,j,1) == 1) then
          ingb = 0
          do m = max(this%is, i - 1), min(this%ie, i + 1)
            if (this%cell_mask(m,1,1) == 1 .and. m /= i) then
              ingb = ingb + 1
              this%cell_ngb(1,ingb,i,j,k) = m
            end if
          end do
          this%cell_ngb_size(i,j,k) = ingb
        else
          this%cell_ngb_size(i,j,k) = 0
        end if
      end do
    end if

    allocate(this%xc(this%is:this%ie))
    if (present(xc)) then
      this%xc = xc
    else
      ! Set coordinates of cells on the large stencil with origin at center.
      do i = this%is, this%ie
        this%xc(i) = -int(sw / 2) + i - 1
      end do
    end if

    ! Initialize sub-stencils.
    if (this%id == 0) then
      this%sub_sw = int((sw + 1) / 2) ! Set sub-stencil size.
      this%sub_nc = this%sub_sw**nd
      this%ns     = this%sub_nc
      this%nterm  = this%ns
      allocate(this%subs(this%ns))
      isub = 1
      j = 1
      do i = this%is, this%is + this%sub_sw - 1
        ! init sub procedures
        this%subs(isub)%init                   => weno_1d_init
        this%subs(isub)%calc_recon_matrix      => weno_1d_calc_recon_matrix     
        this%subs(isub)%calc_ideal_coefs       => weno_1d_calc_ideal_coefs      
        this%subs(isub)%reconstruct            => weno_1d_reconstruct           
        this%subs(isub)%trouble_cell_indicator => weno_1d_trouble_cell_indicator
      
        sub_ie = i + this%sub_sw - 1
        sub_je = 1
        call this%subs(isub)%init(nd=nd, sw=this%sub_sw, xc=this%xc(i:sub_ie), &
                                  dx=this%dx, is=i, ie=sub_ie, id=isub, &
                                  mask=this%cell_mask(i:sub_ie,:,:))
        
        select case (this%sub_sw)
        case (2)
          this%subs(isub)%smooth_indicator => smooth_indicator_2
        case (3)
          this%subs(isub)%smooth_indicator => smooth_indicator_3
        case default 
          stop 'Unknown stencil width in WENO1D, choose only recon_v_order=3 or recon_v_order=5'
        end select
        isub = isub + 1
      end do
      
      select case(this%sw)
      case(3)
        this%get_jab => WENO3_jab
      case(5)
        this%get_jab => WENO5_Z_jab
      case default
        stop 'Unsupported WENO stencil width for get Jacobian, choose stencil width from 3 or 5'
      end select
    end if

    this%initialized = .true.

  end subroutine weno_1d_init

  subroutine weno_1d_calc_recon_matrix(this, ierr)

    class(recon_type), intent(inout) :: this
    integer, intent(out) :: ierr

    ! Local double double arrays for preserving precision.
    real(16), allocatable, dimension(:,:) :: A, iA
    integer , allocatable :: avaCellIdx(:)
    integer i, j, k, ipt, n, iterm, dxn, dyn
    integer iCell, iAvaCell
    
    ierr = 0
    
    if (allocated(this%poly         )) deallocate(this%poly         )
    if (allocated(this%poly_r16     )) deallocate(this%poly_r16     )
    if (allocated(this%iA           )) deallocate(this%iA           )
    if (allocated(this%recon_mtx    )) deallocate(this%recon_mtx    )
    if (allocated(this%recon_mtx_r16)) deallocate(this%recon_mtx_r16)
    if (allocated(this%dpoly_r16    )) deallocate(this%dpoly_r16    )
    if (allocated(this%drecon_mtx   )) deallocate(this%drecon_mtx   )

    allocate(this%poly         (this%npt,this%nc)); this%poly          = 0
    allocate(this%poly_r16     (this%npt,this%nc)); this%poly_r16      = 0
    allocate(this%iA           (this%nc ,this%nc)); this%iA            = 0
    allocate(this%recon_mtx    (this%npt,this%nc)); this%recon_mtx     = 0
    allocate(this%recon_mtx_r16(this%npt,this%nc)); this%recon_mtx_r16 = 0
    allocate(this%dpoly_r16    (this%npt,this%nc,0:this%swx-1,0:this%swy-1,0:this%swz-1))
    allocate(this%drecon_mtx   (this%npt,this%nc,0:this%swx-1,0:this%swy-1,0:this%swz-1))

    allocate( A        (this%nc,this%nc))
    allocate(iA        (this%nc,this%nc)); iA = 0
    allocate(avaCellIdx(this%nc        )); avaCellIdx = 0

    ! Set index maps.
    iAvaCell = 1
    iCell = 1
    j = 1
    do i = this%is, this%ie
      if (this%cell_mask(i,j,1) == 1) then
        avaCellIdx(iCell) = iAvaCell
        iAvaCell = iAvaCell + 1
      end if
      iCell = iCell + 1
    end do

    ! Set the poly for each evaluation point.
    ! Select monomials according to mask.
    do ipt = 1, this%npt
      iterm = 1
      do i = 1, this%sw
        if (this%poly_mask(i,1,1) == 1) then
          call calc_monomial(this%x(ipt), i - 1, this%poly_r16(ipt,iterm))
          iterm = iterm + 1
        end if
      end do
    end do
    
    do ipt = 1, this%npt
      iterm = 1
      do i = 1, this%sw
        do dxn = 0, this%sw - 1
          call calc_deriv_monomial(this%x(ipt), iterm - 1, 1, this%sw, this%dpoly_r16(ipt,iterm,dxn,0,0))
          this%dpoly_r16(ipt,iterm,dxn,0,0) = this%dpoly_r16(ipt,iterm,dxn,0,0) / this%dx**dxn
        end do
        iterm = iterm + 1
      end do
    end do

    ! Calculate inverse of integral coefficient matrix.
    call calc_tensor_product_poly_integral_matrix(this%sw, this%xc, A, this%cell_mask(:,1,1), this%poly_mask(:,1,1))
    n = count(this%cell_mask == 1)
    call inverse_matrix(A(1:n,1:n), iA(1:n,1:n), ierr)
    if (ierr /= 0) then
      print*,'Failed to calculate inverse A in weno_1d_calc_recon_matrix'
      deallocate(A, iA, avaCellIdx)
      return
    end if

    this%recon_mtx_r16(:,1:n) = matmul(this%poly_r16(:,1:n), iA(1:n,1:n))

    ! Copy double double iA into double iA.
    this%iA   = iA
    this%poly = this%poly_r16
    
    do dxn = 0, this%swx - 1
      this%drecon_mtx(1:this%npt,1:n,dxn,0,0) = matmul(this%dpoly_r16(1:this%npt,1:n,dxn,0,0), iA(1:n,1:n))
    end do

    ! Rearrange iA and recon_mtx.
    do n = this%nc, 1, -1
      if (avaCellIdx(n) /= 0) then
        this%iA           (:,n) = this%iA           (:,avaCellIdx(n))
        this%recon_mtx_r16(:,n) = this%recon_mtx_r16(:,avaCellIdx(n))
      else
        this%iA           (:,n) = 0
        this%recon_mtx_r16(:,n) = 0
      end if
    end do

    ! Copy double double recon_mtx into double recon_mtx.
    this%recon_mtx = this%recon_mtx_r16

    deallocate(A, iA, avaCellIdx)

  end subroutine weno_1d_calc_recon_matrix

  subroutine weno_1d_calc_ideal_coefs(this, ierr)

    class(recon_type), intent(inout) :: this
    integer, intent(out) :: ierr

    ! Local double double arrays for preserving precision.
    real(16), allocatable, dimension(:,:) :: A, ATA, iATA, gamma
    integer i, j, k, ipt, m, n, di, dj, isub

    ierr = 0

    if (.not. this%initialized) then
      ierr = 1
      return
    end if
    
    do isub = 1, this%ns
      call this%subs(isub)%calc_recon_matrix(ierr)
      if (ierr /= 0) return
    end do
    call this%calc_recon_matrix(ierr)
    if (ierr /= 0) return

    if (allocated(this%gamma)) deallocate(this%gamma)
    allocate(this%gamma(this%ns,this%npt))

    allocate(    A(this%nc,this%ns )); A = 0
    allocate(  ATA(this%ns,this%ns ))
    allocate( iATA(this%ns,this%ns ))
    allocate(gamma(this%ns,this%npt)); gamma = 0
    do ipt = 1, this%npt
      do isub = 1, this%ns
        di = this%subs(isub)%is - this%is
        dj = 0
        j = 1
        do i = 1, this%sub_sw
          if (this%subs(isub)%cell_mask(this%subs(isub)%is+i-1,1,1) == 1) then
            m = (j + dj - 1) * this%sw     + i + di
            n = (j      - 1) * this%sub_sw + i
            A(m,isub) = this%subs(isub)%recon_mtx_r16(ipt,n)
          end if
        end do
      end do
      ATA = matmul(transpose(A), A)
      call inverse_matrix(ATA, iATA, ierr)
      if (ierr /= 0) return
      gamma(:,ipt) = matmul(matmul(iATA, transpose(A)), this%recon_mtx_r16(ipt,:))
      ! Set near-zero values to zero exactly.
      do isub = 1, this%ns
        if (abs(gamma(isub,ipt)) < 1.0e-15) gamma(isub,ipt) = 0
      end do
    end do
    if (abs(sum(matmul(A, gamma(:,1)) - this%recon_mtx_r16(1,:))) > 1.0e-15) then
      ierr = 3
      this%gamma = 0
    else
      this%gamma = gamma
    end if
    deallocate(A, ATA, iATA, gamma)

  end subroutine weno_1d_calc_ideal_coefs

  subroutine weno_1d_reconstruct(this, fi, fo, TCI, ierr)

    class(recon_type), intent(inout) :: this
    real(8), intent(in ) :: fi(:,:,:) ! Cell averaged function values
    real(8), intent(out) :: fo(:)     ! Reconstructed function values on evaluation points
    logical, intent(in ), optional :: TCI ! Trouble Cell Indicator(TCI), 1 for existing TC(use WENO), 0 for no TC(use poly)
    integer, intent(out), optional :: ierr

    integer isub, ipt
    real(8) a     (this%nterm,this%ns)
    real(8) beta  (this%ns)
    real(8) alpha (this%ns)
    real(8) w     (this%ns)
    
    real(8) fs(this%ns,this%npt)
    
    ! Calculate point values from each available sub-stencils and smoothness
    ! indicators for those sub-stencils.
    do isub = 1, this%ns
      a(:,isub) = matmul(this%subs(isub)%iA, pack(fi(this%subs(isub)%is:this%subs(isub)%ie,1,1), .true.))
      beta(isub) = this%subs(isub)%smooth_indicator(a(:,isub))
      fs(isub,:) = matmul(this%subs(isub)%poly, a(:,isub))
    end do

    alpha = weno_z_alpha(beta, this%ns, this%eps)

    do ipt = 1, this%npt
      w = this%gamma(:,ipt) * alpha
      fo(ipt) = sum( w * fs(:,ipt) ) / sum(w)
    end do
    
  end subroutine weno_1d_reconstruct

  pure function weno_z_alpha(beta, ns, eps) result(res)

    real(8), intent(in) :: beta(ns)
    integer, intent(in) :: ns
    real(8), intent(in) :: eps
    real(8) res(ns)

    real(8) tau
    real(8) xi(ns)

    tau = abs( beta(ns) - beta(1) )
    xi = ( tau + eps ) / ( beta + eps )
    
    res = 1. + xi**2

  end function weno_z_alpha

  pure logical function weno_1d_trouble_cell_indicator(this, fi, dx) result(res)

    class(recon_type), intent(in) :: this
    real(8), intent(in) :: fi(:,:,:)
    real(8), intent(in) :: dx

    integer i, j, k, m, n, ingb
    real(8), parameter :: sqrt2 = sqrt(2.0d0)
    real(8) c
    
    res = .true.
    !stop 'No TCI support in WENO1D' 
    
  end function weno_1d_trouble_cell_indicator
  
  subroutine weno_1d_reconstruct_deriv(this, fi, fx, fy, fz, ierr)
    class(recon_type), intent(inout) :: this
    real(8), intent(in ) :: fi(:,:,:) ! Cell averaged function values
    real(8), intent(out), optional :: fx(:)     ! Reconstructed x derivative values on evaluation points
    real(8), intent(out), optional :: fy(:)     ! Reconstructed y derivative values on evaluation points
    real(8), intent(out), optional :: fz(:)     ! Reconstructed y derivative values on evaluation points
    integer, intent(out), optional :: ierr

#ifdef DEBUG
    if(present(ierr))then
      ierr = 0
      if (size(fi) /= this%nc) then
        if (present(ierr)) ierr = 1
        return
      end if
      if (size(fo) /= this%npt) then
        if (present(ierr)) ierr = 1
        return
      end if
    endif
#endif

    fx = matmul(this%drecon_mtx(:,:,1,0,0), pack(fi, .true.))
    
  end subroutine weno_1d_reconstruct_deriv
  
  subroutine WENO3_jab(this,dqLdq,dqRdq,q)
    class(recon_type), intent(in) :: this
    real(8), dimension(this%sw), intent(inout) :: dqLdq,dqRdq
    real(8), dimension(this%sw), intent(in   ), optional :: q
    
    real(8) :: q1,q2,q3,q4,q5
    
    real(8), parameter :: c1 = 1./3.
    real(8), parameter :: c2 = 2./3.
    
    real(8) :: &
    t2   ,&
    t3   ,&
    t4   ,&
    t5   ,&
    t7   ,&
    t9   ,&
    t10  ,&
    t11  ,&
    t12  ,&
    t6   ,&
    t8   ,&
    t13  ,&
    t14  ,&
    t15  ,&
    t16  ,&
    t24  ,&
    t25  ,&
    t17  ,&
    t18  ,&
    t19  ,&
    t20  ,&
    t21  ,&
    t27  ,&
    t28  ,&
    t22  ,&
    t23  ,&
    t26  ,&
    t29  ,&
    t30  ,&
    t31  ,&
    t32  ,&
    t33  ,&
    t34  ,&
    t35  ,&
    t36  ,&
    t47  ,&
    t48  ,&
    t49  ,&
    t50  ,&
    t51  ,&
    t52  ,&
    t53  ,&
    t54  ,&
    t37  ,&
    t38  ,&
    t45  ,&
    t46  ,&
    t55  ,&
    t56  ,&
    t57  ,&
    t58  ,&
    t39  ,&
    t40  ,&
    t65  ,&
    t66  ,&
    t67  ,&
    t68  ,&
    t41  ,&
    t42  ,&
    t43  ,&
    t44  ,&
    t69  ,&
    t70  ,&
    t71  ,&
    t72  ,&
    t73  ,&
    t74  ,&
    t75  ,&
    t76  ,&
    t59  ,&
    t60  ,&
    t77  ,&
    t78  ,&
    t81  ,&
    t82  ,&
    t83  ,&
    t84  ,&
    t61  ,&
    t62  ,&
    t85  ,&
    t86  ,&
    t63  ,&
    t64  ,&
    t79  ,&
    t80
    
    q1 = q(1)
    q2 = q(2)
    q3 = q(3)
    
    t2 = q1*2.0;
    t3 = q2*2.0;
    t4 = q3*2.0;
    t5 = -q2;
    t7 = -q3;
    t9 = q1/2.0;
    t10 = q2/2.0;
    t11 = q3/2.0;
    t12 = q2*(3.0/2.0);
    t6 = -t3;
    t8 = -t4;
    t13 = q1+t5;
    t14 = q2+t7;
    t15 = -t11;
    t16 = -t12;
    t24 = t9+t10;
    t25 = t10+t11;
    t17 = t2+t6;
    t18 = t2+t8;
    t19 = t3+t8;
    t20 = t13*t13;
    t21 = t14*t14;
    t27 = t9+t16;
    t28 = t12+t15;
    t22 = this%eps+t20;
    t23 = this%eps+t21;
    t26 = -t21;
    t29 = 1.0/(t22*t22);
    t30 = 1.0/(t22*t22*t22);
    t31 = 1.0/(t23*t23);
    t32 = 1.0/(t23*t23*t23);
    t33 = t20+t26;
    t34 = abs(t33);
    t35 = sign(1._8,t33)!(t33/abs(t33));
    t36 = t34*t34;
    t47 = t17*t29*t34*t35*2.0;
    t48 = t18*t29*t34*t35*2.0;
    t49 = t18*t31*t34*t35*2.0;
    t50 = t19*t31*t34*t35*2.0;
    t51 = C1*t17*t31*t34*t35*2.0;
    t52 = C1*t19*t29*t34*t35*2.0;
    t53 = C2*t17*t31*t34*t35*2.0;
    t54 = C2*t19*t29*t34*t35*2.0;
    t37 = t29*t36;
    t38 = t31*t36;
    t45 = t17*t30*t36*2.0;
    t46 = t19*t32*t36*2.0;
    t55 = -t47;
    t56 = -t48;
    t57 = -t51;
    t58 = -t53;
    t39 = t37+1.0;
    t40 = t38+1.0;
    t65 = t46+t49;
    t66 = t46+t50;
    t67 = t45+t55;
    t68 = t45+t56;
    t41 = C1*t39;
    t42 = C2*t39;
    t43 = C1*t40;
    t44 = C2*t40;
    t69 = C1*t65;
    t70 = C1*t66;
    t71 = C2*t65;
    t72 = C2*t66;
    t73 = C1*t67;
    t74 = C1*t68;
    t75 = C2*t67;
    t76 = C2*t68;
    t59 = t41+t44;
    t60 = t42+t43;
    t77 = -t69;
    t78 = -t71;
    t81 = t54+t70;
    t82 = t52+t72;
    t83 = t58+t73;
    t84 = t57+t75;
    t61 = 1.0/t59;
    t62 = 1.0/t60;
    t85 = t74+t78;
    t86 = t76+t77;
    t63 = t61*t61;
    t64 = t62*t62;
    t79 = (t42*t62)/2.0;
    t80 = (t44*t61)/2.0;
    
    dqLdq(1) = t79-t24*t62*t75+t43*t64*(t11+t16)*(t51-t75)-t24*t42*t64*(t51-t75)-C1*t17*t31*t34*t35*t62*(t11+t16)*2.0;
    dqLdq(2) = t79+t43*t62*(3.0/2.0)+t62*t69*(t11+t16)+t24*t62*t76-t43*t64*(t11+t16)*(t69-t76)+t24*t42*t64*(t69-t76);
    dqLdq(3) = t43*t62*(-1.0/2.0)-t62*t70*(t11+t16)+t24*t54*t62+t43*t64*t81*(t11+t16)-t24*t42*t64*t81;
    dqRdq(1) = t41*t61*(-1.0/2.0)+t25*t53*t61+t27*t61*t73+t27*t41*t63*(t53-t73)-t25*t44*t63*(t53-t73);
    dqRdq(2) = t80+t41*t61*(3.0/2.0)-t27*t61*t74+t25*t61*t78-t27*t41*t63*(t71-t74)+t25*t44*t63*(t71-t74);
    dqRdq(3) = t80+t25*t61*t72+t27*t41*t63*t82-t25*t44*t63*t82-C1*t19*t27*t29*t34*t35*t61*2.0;
    
  end subroutine WENO3_jab
  
  subroutine WENO5_Z_jab(this,dqLdq,dqRdq,q)
    class(recon_type), intent(in) :: this
    real(8), dimension(this%nc), intent(inout) :: dqLdq,dqRdq
    real(8), dimension(this%nc), intent(in   ), optional :: q
    
    real(8) :: q1,q2,q3,q4,q5
    
    real(8), parameter :: c1      = 0.1
    real(8), parameter :: c2      = 0.6
    real(8), parameter :: c3      = 0.3
    
    real(8) :: &
    t2  ,&
    t3  ,&
    t4  ,&
    t5  ,&
    t6  ,&
    t7  ,&
    t10  ,&
    t14  ,&
    t15  ,&
    t16  ,&
    t17  ,&
    t18  ,&
    t19  ,&
    t20  ,&
    t21  ,&
    t22  ,&
    t23  ,&
    t24  ,&
    t25  ,&
    t26  ,&
    t27  ,&
    t28  ,&
    t29  ,&
    t34  ,&
    t36  ,&
    t37  ,&
    t39  ,&
    t40  ,&
    t41  ,&
    t43  ,&
    t44  ,&
    t45  ,&
    t46  ,&
    t47  ,&
    t48  ,&
    t49  ,&
    t50  ,&
    t51  ,&
    t52  ,&
    t53  ,&
    t54  ,&
    t55  ,&
    t56  ,&
    t57  ,&
    t58  ,&
    t8  ,&
    t9  ,&
    t11  ,&
    t12  ,&
    t13  ,&
    t30  ,&
    t31  ,&
    t32  ,&
    t33  ,&
    t35  ,&
    t38  ,&
    t42  ,&
    t62  ,&
    t63  ,&
    t64  ,&
    t65  ,&
    t66  ,&
    t67  ,&
    t68  ,&
    t69  ,&
    t70  ,&
    t59  ,&
    t60  ,&
    t61  ,&
    t71  ,&
    t72  ,&
    t73  ,&
    t87  ,&
    t88  ,&
    t89  ,&
    t90  ,&
    t91  ,&
    t92  ,&
    t93  ,&
    t94  ,&
    t95  ,&
    t96  ,&
    t97  ,&
    t98  ,&
    t99  ,&
    t100 ,&
    t101 ,&
    t102 ,&
    t74  ,&
    t75  ,&
    t76  ,&
    t77  ,&
    t78  ,&
    t79  ,&
    t80  ,&
    t81  ,&
    t82  ,&
    t83  ,&
    t84  ,&
    t85  ,&
    t86  ,&
    t103 ,&
    t106 ,&
    t107 ,&
    t104 ,&
    t105 ,&
    t108 ,&
    t109 ,&
    t110 ,&
    t111 ,&
    t112 ,&
    t113 ,&
    t114 ,&
    t115 ,&
    t116 ,&
    t151 ,&
    t152 ,&
    t153 ,&
    t154 ,&
    t156 ,&
    t157 ,&
    t158 ,&
    t159 ,&
    t160 ,&
    t161 ,&
    t162 ,&
    t163 ,&
    t164 ,&
    t165 ,&
    t166 ,&
    t167 ,&
    t168 ,&
    t169 ,&
    t170 ,&
    t171 ,&
    t172 ,&
    t173 ,&
    t174 ,&
    t175 ,&
    t176 ,&
    t117 ,&
    t118 ,&
    t119 ,&
    t120 ,&
    t121 ,&
    t127 ,&
    t128 ,&
    t130 ,&
    t131 ,&
    t132 ,&
    t133 ,&
    t134 ,&
    t135 ,&
    t136 ,&
    t137 ,&
    t138 ,&
    t139 ,&
    t140 ,&
    t141 ,&
    t142 ,&
    t155 ,&
    t122 ,&
    t123 ,&
    t124 ,&
    t125 ,&
    t126 ,&
    t129 ,&
    t143 ,&
    t144 ,&
    t145 ,&
    t146 ,&
    t147 ,&
    t148 ,&
    t149 ,&
    t150 ,&
    t177 ,&
    t178 ,&
    t179 ,&
    t180 ,&
    t181 ,&
    t186 ,&
    t187 ,&
    t182 ,&
    t183 ,&
    t188 ,&
    t189 ,&
    t190 ,&
    t191 ,&
    t192 ,&
    t193 ,&
    t194 ,&
    t195 ,&
    t184 ,&
    t185
    
    q1 = q(1)
    q2 = q(2)
    q3 = q(3)
    q4 = q(4)
    q5 = q(5)
    
    t2 = q2*2.0
    t3 = q3*2.0
    t4 = q2*4.0
    t5 = q3*3.0
    t6 = q4*2.0
    t7 = q4*4.0
    t10 = -q4
    t14 = q1/3.0
    t15 = q2/3.0
    t16 = q3/3.0
    t17 = q1/6.0
    t18 = q4/3.0
    t19 = q2/6.0
    t20 = q5/3.0
    t21 = q2*(5.0/3.0)
    t22 = q4/6.0
    t23 = q4*(5.0/3.0)
    t24 = q5/6.0
    t25 = q2*(5.0/6.0)
    t26 = q3*(5.0/6.0)
    t27 = q2*(7.0/6.0)
    t28 = q4*(5.0/6.0)
    t29 = q4*(7.0/6.0)
    t34 = q1*(8.0/3.0)
    t36 = q2*(8.0/3.0)
    t37 = q1*(1.1E+1/3.0)
    t39 = q4*(8.0/3.0)
    t40 = q5*(8.0/3.0)
    t41 = q3*(1.1E+1/3.0)
    t43 = q2*(1.3E+1/3.0)
    t44 = q3*(1.3E+1/3.0)
    t45 = q5*(1.1E+1/3.0)
    t46 = q3*(1.1E+1/6.0)
    t47 = q4*(1.3E+1/3.0)
    t48 = q1*(1.9E+1/3.0)
    t49 = q2*(1.9E+1/3.0)
    t50 = q3*(2.0E+1/3.0)
    t51 = q4*(1.9E+1/3.0)
    t52 = q5*(1.9E+1/3.0)
    t53 = q3*(2.6E+1/3.0)
    t54 = q2*(3.1E+1/3.0)
    t55 = q3*(3.1E+1/3.0)
    t56 = q4*(3.1E+1/3.0)
    t57 = q2*(5.0E+1/3.0)
    t58 = q4*(5.0E+1/3.0)
    t8 = -t2
    t9 = -t3
    t11 = -t4
    t12 = -t6
    t13 = -t7
    t30 = q2+t10
    t31 = -t17
    t32 = -t19
    t33 = -t22
    t35 = -t24
    t38 = -t27
    t42 = -t29
    t62 = -t44
    t63 = -t45
    t64 = -t49
    t65 = -t51
    t66 = -t53
    t67 = -t54
    t68 = -t56
    t69 = -t57
    t70 = -t58
    t59 = q1+q3+t8
    t60 = q2+q4+t9
    t61 = q3+q5+t12
    t71 = t30*t30
    t72 = q1+t5+t11
    t73 = q5+t5+t13
    t87 = t16+t25+t31
    t88 = t15+t26+t33
    t89 = t18+t26+t32
    t90 = t16+t28+t35
    t91 = t14+t38+t46
    t92 = t20+t42+t46
    t93 = t21+t39+t62
    t94 = t23+t36+t62
    t95 = t34+t41+t64
    t96 = t40+t41+t65
    t97 = t43+t47+t66
    t98 = t37+t50+t67
    t99 = t45+t50+t68
    t100 = t48+t55+t69
    t101 = t52+t55+t70
    t102 = t37+t56+t63+t67
    t74 = t59*t59
    t75 = t60*t60
    t76 = t61*t61
    t77 = t72*t72
    t78 = t73*t73
    t79 = t71/4.0
    t80 = t77/4.0
    t81 = t78/4.0
    t82 = t74*(1.3E+1/1.2E+1)
    t83 = t75*(1.3E+1/1.2E+1)
    t84 = t76*(1.3E+1/1.2E+1)
    t85 = -t81
    t86 = -t84
    t103 = this%eps+t79+t83
    t106 = this%eps+t80+t82
    t107 = this%eps+t81+t84
    t104 = 1.0/(t103*t103)
    t105 = 1.0/(t103*t103*t103)
    t108 = 1.0/(t106*t106)
    t109 = 1.0/(t106*t106*t106)
    t110 = 1.0/(t107*t107)
    t111 = 1.0/(t107*t107*t107)
    t112 = t80+t82+t85+t86
    t113 = abs(t112)
    t114 = merge( 0._8, sign(1._8,t112), t112==0 )!(t112/abs(t112))
    t115 = this%eps+t113
    t116 = t115*t115
    t151 = t95*t104*t114*t115*(6.0/5.0)
    t152 = t96*t104*t114*t115*(6.0/5.0)
    t153 = t100*t104*t114*t115*(6.0/5.0)
    t154 = t101*t104*t114*t115*(6.0/5.0)
    t156 = (t95*t108*t114*t115)/5.0
    t157 = t95*t108*t114*t115*(3.0/5.0)
    t158 = (t96*t108*t114*t115)/5.0
    t159 = t96*t108*t114*t115*(3.0/5.0)
    t160 = (t95*t110*t114*t115)/5.0
    t161 = t95*t110*t114*t115*(3.0/5.0)
    t162 = (t96*t110*t114*t115)/5.0
    t163 = t96*t110*t114*t115*(3.0/5.0)
    t164 = (t100*t108*t114*t115)/5.0
    t165 = t100*t108*t114*t115*(3.0/5.0)
    t166 = (t101*t108*t114*t115)/5.0
    t167 = t101*t108*t114*t115*(3.0/5.0)
    t168 = (t100*t110*t114*t115)/5.0
    t169 = t100*t110*t114*t115*(3.0/5.0)
    t170 = (t101*t110*t114*t115)/5.0
    t171 = t101*t110*t114*t115*(3.0/5.0)
    t172 = t102*t104*t114*t115*(6.0/5.0)
    t173 = (t102*t108*t114*t115)/5.0
    t174 = t102*t108*t114*t115*(3.0/5.0)
    t175 = (t102*t110*t114*t115)/5.0
    t176 = t102*t110*t114*t115*(3.0/5.0)
    t117 = t104*t116*(3.0/5.0)
    t118 = (t108*t116)/1.0E+1
    t119 = t108*t116*(3.0/1.0E+1)
    t120 = (t110*t116)/1.0E+1
    t121 = t110*t116*(3.0/1.0E+1)
    t127 = t93*t105*t116*(6.0/5.0)
    t128 = t94*t105*t116*(6.0/5.0)
    t130 = t97*t105*t116*(6.0/5.0)
    t131 = (t95*t109*t116)/5.0
    t132 = t95*t109*t116*(3.0/5.0)
    t133 = (t96*t111*t116)/5.0
    t134 = t96*t111*t116*(3.0/5.0)
    t135 = (t98*t109*t116)/5.0
    t136 = t98*t109*t116*(3.0/5.0)
    t137 = (t99*t111*t116)/5.0
    t138 = t99*t111*t116*(3.0/5.0)
    t139 = (t100*t109*t116)/5.0
    t140 = t100*t109*t116*(3.0/5.0)
    t141 = (t101*t111*t116)/5.0
    t142 = t101*t111*t116*(3.0/5.0)
    t155 = -t154
    t122 = t117+3.0/5.0
    t123 = t118+1.0/1.0E+1
    t124 = t119+3.0/1.0E+1
    t125 = t120+1.0/1.0E+1
    t126 = t121+3.0/1.0E+1
    t129 = -t127
    t143 = -t131
    t144 = -t132
    t145 = -t135
    t146 = -t136
    t147 = -t137
    t148 = -t138
    t149 = -t139
    t150 = -t140
    t177 = t128+t153
    t178 = t127+t155
    t179 = t130+t172
    t180 = t117+t118+t121+1.0
    t181 = t117+t119+t120+1.0
    t186 = t133+t152+t159+t162
    t187 = t134+t152+t158+t163
    t182 = 1.0/t180
    t183 = 1.0/t181
    t188 = t143+t151+t156+t161
    t189 = t144+t151+t157+t160
    t190 = t149+t164+t169+t177
    t191 = t150+t165+t168+t177
    t192 = t129+t141+t154+t167+t170
    t193 = t129+t142+t154+t166+t171
    t194 = t145+t148+t173+t176+t179
    t195 = t146+t147+t174+t175+t179
    t184 = t182*t182
    t185 = t183*t183
    
    dqLdq(1) = t124*t183*(-1.0/6.0)+t88*t151*t183+t92*t160*t183-t87*t183*(t132-t157)-t88*t122*t185*t189-t87*t124*t185*t189-t92*t125*t185*t189
    dqLdq(2) = (t122*t183)/3.0+t124*t183*(5.0/6.0)-t88*t177*t183+t87*t183*(t140-t165)+t88*t122*t185*t191+t87*t124*t185*t191+t92*t125*t185*t191-(t92*t100*t110*t114*t115*t183)/5.0
    dqLdq(3) = t122*t183*(5.0/6.0)+(t124*t183)/3.0+t125*t183*(1.1E+1/6.0)+t88*t179*t183-t87*t183*(t136-t174)-t92*t183*(t137-t175)-t88*t122*t185*t195-t87*t124*t185*t195-t92*t125*t185*t195
    dqLdq(4) = t122*t183*(-1.0/6.0)-t125*t183*(7.0/6.0)+t92*t183*(t141+t170)+t87*t167*t183-t88*t178*t183-t88*t122*t185*t192-t87*t124*t185*t192-t92*t125*t185*t192
    dqLdq(5) = (t125*t183)/3.0-t92*t183*(t133+t162)+t88*t122*t185*t186+t87*t124*t185*t186+t92*t125*t185*t186-t88*t96*t104*t114*t115*t183*(6.0/5.0)-t87*t96*t108*t114*t115*t183*(3.0/5.0)
    dqRdq(1) = (t123*t182)/3.0+t89*t151*t182+t90*t161*t182-t91*t182*(t131-t156)-t89*t122*t184*t188-t91*t123*t184*t188-t90*t126*t184*t188
    dqRdq(2) = t122*t182*(-1.0/6.0)-t123*t182*(7.0/6.0)-t89*t177*t182+t91*t182*(t139-t164)+t89*t122*t184*t190+t91*t123*t184*t190+t90*t126*t184*t190-t90*t100*t110*t114*t115*t182*(3.0/5.0)
    dqRdq(3) = t122*t182*(5.0/6.0)+t123*t182*(1.1E+1/6.0)+(t126*t182)/3.0+t89*t179*t182-t91*t182*(t135-t173)-t90*t182*(t138-t176)-t89*t122*t184*t194-t91*t123*t184*t194-t90*t126*t184*t194
    dqRdq(4) = (t122*t182)/3.0+t126*t182*(5.0/6.0)+t90*t182*(t142+t171)+t91*t166*t182-t89*t178*t182-t89*t122*t184*t193-t91*t123*t184*t193-t90*t126*t184*t193
    dqRdq(5) = t126*t182*(-1.0/6.0)-t90*t182*(t134+t163)+t89*t122*t184*t187+t91*t123*t184*t187+t90*t126*t184*t187-t89*t96*t104*t114*t115*t182*(6.0/5.0)-(t91*t96*t108*t114*t115*t182)/5.0
      
  end subroutine WENO5_Z_jab
  
  subroutine WENO5_JS_jab(this,dqLdq,dqRdq,q)
    class(recon_type), intent(in) :: this
    real(8), dimension(this%nc), intent(inout) :: dqLdq,dqRdq
    real(8), dimension(this%nc), intent(in   ) :: q
    
    real(8) :: q1,q2,q3,q4,q5
    
    real(8), parameter :: c1      = 0.1
    real(8), parameter :: c2      = 0.6
    real(8), parameter :: c3      = 0.3
    
    real(8) :: &
    t2     ,&
    t3     ,&
    t4     ,&
    t5     ,&
    t6     ,&
    t7     ,&
    t10     ,&
    t14     ,&
    t15     ,&
    t16     ,&
    t17     ,&
    t18     ,&
    t19     ,&
    t20     ,&
    t21     ,&
    t22     ,&
    t23     ,&
    t24     ,&
    t25     ,&
    t26     ,&
    t27     ,&
    t28     ,&
    t29     ,&
    t34     ,&
    t36     ,&
    t37     ,&
    t39     ,&
    t40     ,&
    t41     ,&
    t43     ,&
    t44     ,&
    t45     ,&
    t46     ,&
    t47     ,&
    t48     ,&
    t49     ,&
    t50     ,&
    t51     ,&
    t52     ,&
    t53     ,&
    t54     ,&
    t55     ,&
    t56     ,&
    t57     ,&
    t58     ,&
    t8     ,&
    t9     ,&
    t11     ,&
    t12     ,&
    t13     ,&
    t30     ,&
    t31     ,&
    t32     ,&
    t33     ,&
    t35     ,&
    t38     ,&
    t42     ,&
    t62     ,&
    t63     ,&
    t64     ,&
    t65     ,&
    t66     ,&
    t67     ,&
    t68     ,&
    t69     ,&
    t59     ,&
    t60     ,&
    t61     ,&
    t70     ,&
    t71     ,&
    t72     ,&
    t84     ,&
    t85     ,&
    t86     ,&
    t87     ,&
    t88     ,&
    t89     ,&
    t90     ,&
    t91     ,&
    t92     ,&
    t93     ,&
    t94     ,&
    t95     ,&
    t96     ,&
    t97     ,&
    t98     ,&
    t73     ,&
    t74     ,&
    t75     ,&
    t76     ,&
    t77     ,&
    t78     ,&
    t79     ,&
    t80     ,&
    t81     ,&
    t82     ,&
    t83     ,&
    t99     ,&
    t102    ,&
    t103    ,&
    t100    ,&
    t101    ,&
    t104    ,&
    t105    ,&
    t106    ,&
    t107    ,&
    t108    ,&
    t109    ,&
    t110    ,&
    t111    ,&
    t112    ,&
    t113    ,&
    t114    ,&
    t115    ,&
    t116    ,&
    t119    ,&
    t121    ,&
    t122    ,&
    t123    ,&
    t124    ,&
    t125    ,&
    t126    ,&
    t127    ,&
    t128    ,&
    t117    ,&
    t118    ,&
    t120    ,&
    t133    ,&
    t134    ,&
    t129    ,&
    t130    ,&
    t131    ,&
    t132    ,&
    t135    ,&
    t136    ,&
    t139    ,&
    t140    ,&
    t137    ,&
    t138
    
    q1 = q(1)
    q2 = q(2)
    q3 = q(3)
    q4 = q(4)
    q5 = q(5)
    
    t2 = q2*2.0
    t3 = q3*2.0
    t4 = q2*4.0
    t5 = q3*3.0
    t6 = q4*2.0
    t7 = q4*4.0
    t10 = -q4
    t14 = q1/3.0
    t15 = q2/3.0
    t16 = q3/3.0
    t17 = q1/6.0
    t18 = q4/3.0
    t19 = q2/6.0
    t20 = q5/3.0
    t21 = q2*(5.0/3.0)
    t22 = q4/6.0
    t23 = q4*(5.0/3.0)
    t24 = q5/6.0
    t25 = q2*(5.0/6.0)
    t26 = q3*(5.0/6.0)
    t27 = q2*(7.0/6.0)
    t28 = q4*(5.0/6.0)
    t29 = q4*(7.0/6.0)
    t34 = q1*(8.0/3.0)
    t36 = q2*(8.0/3.0)
    t37 = q1*(1.1E+1/3.0)
    t39 = q4*(8.0/3.0)
    t40 = q5*(8.0/3.0)
    t41 = q3*(1.1E+1/3.0)
    t43 = q2*(1.3E+1/3.0)
    t44 = q3*(1.3E+1/3.0)
    t45 = q5*(1.1E+1/3.0)
    t46 = q3*(1.1E+1/6.0)
    t47 = q4*(1.3E+1/3.0)
    t48 = q1*(1.9E+1/3.0)
    t49 = q2*(1.9E+1/3.0)
    t50 = q3*(2.0E+1/3.0)
    t51 = q4*(1.9E+1/3.0)
    t52 = q5*(1.9E+1/3.0)
    t53 = q3*(2.6E+1/3.0)
    t54 = q2*(3.1E+1/3.0)
    t55 = q3*(3.1E+1/3.0)
    t56 = q4*(3.1E+1/3.0)
    t57 = q2*(5.0E+1/3.0)
    t58 = q4*(5.0E+1/3.0)
    t8 = -t2
    t9 = -t3
    t11 = -t4
    t12 = -t6
    t13 = -t7
    t30 = q2+t10
    t31 = -t17
    t32 = -t19
    t33 = -t22
    t35 = -t24
    t38 = -t27
    t42 = -t29
    t62 = -t44
    t63 = -t49
    t64 = -t51
    t65 = -t53
    t66 = -t54
    t67 = -t56
    t68 = -t57
    t69 = -t58
    t59 = q1+q3+t8
    t60 = q2+q4+t9
    t61 = q3+q5+t12
    t70 = t30*t30
    t71 = q1+t5+t11
    t72 = q5+t5+t13
    t84 = t16+t25+t31
    t85 = t15+t26+t33
    t86 = t18+t26+t32
    t87 = t16+t28+t35
    t88 = t14+t38+t46
    t89 = t20+t42+t46
    t90 = t21+t39+t62
    t91 = t23+t36+t62
    t92 = t34+t41+t63
    t93 = t40+t41+t64
    t94 = t43+t47+t65
    t95 = t37+t50+t66
    t96 = t45+t50+t67
    t97 = t48+t55+t68
    t98 = t52+t55+t69
    t73 = t59*t59
    t74 = t60*t60
    t75 = t61*t61
    t76 = t71*t71
    t77 = t72*t72
    t78 = t70/4.0
    t79 = t76/4.0
    t80 = t77/4.0
    t81 = t73*(1.3E+1/1.2E+1)
    t82 = t74*(1.3E+1/1.2E+1)
    t83 = t75*(1.3E+1/1.2E+1)
    t99 = this%eps+t78+t82
    t102 = this%eps+t79+t81
    t103 = this%eps+t80+t83
    t100 = 1.0/(t99*t99)
    t101 = 1.0/(t99*t99*t99)
    t104 = 1.0/(t102*t102)
    t105 = 1.0/(t102*t102*t102)
    t106 = 1.0/(t102*t102*t102*t102*t102)
    t107 = 1.0/(t103*t103)
    t108 = 1.0/(t103*t103*t103)
    t109 = 1.0/(t103*t103*t103*t103*t103)
    t110 = t100*(3.0/5.0)
    t111 = t104/1.0E+1
    t112 = t104*(3.0/1.0E+1)
    t113 = t107/1.0E+1
    t114 = t107*(3.0/1.0E+1)
    t115 = t90*t101*(6.0/5.0)
    t116 = t91*t101*(6.0/5.0)
    t119 = t94*t101*(6.0/5.0)
    t121 = (t95*t105)/5.0
    t122 = t95*t105*(3.0/5.0)
    t123 = (t96*t108)/5.0
    t124 = t96*t108*(3.0/5.0)
    t125 = (t97*t105)/5.0
    t126 = t97*t105*(3.0/5.0)
    t127 = (t98*t108)/5.0
    t128 = t98*t108*(3.0/5.0)
    t117 = -t115
    t118 = -t116
    t120 = -t119
    t133 = t110+t111+t114
    t134 = t110+t112+t113
    t129 = t118+t125
    t130 = t118+t126
    t131 = t117+t127
    t132 = t117+t128
    t135 = 1.0/t133
    t136 = 1.0/t134
    t139 = t120+t121+t124
    t140 = t120+t122+t123
    t137 = t135*t135
    t138 = t136*t136
    
    dqLdq(1) = t104*t136*(-1.0/2.0E+1)-t84*t92*t105*t136*(3.0/5.0)+t84*t92*t106*t138*(9.0/5.0E+1)+t85*t92*t100*t105*t138*(9.0/2.5E+1)+t89*t92*t105*t107*t138*(3.0/5.0E+1)
    dqLdq(2) = (t100*t136)/5.0+(t104*t136)/4.0+t84*t126*t136-t85*t91*t101*t136*(6.0/5.0)+t85*t110*t138*(t116-t126)+t84*t112*t138*(t116-t126)+t89*t113*t138*(t116-t126)
    dqLdq(3) = (t100*t136)/2.0+t107*t136*(1.1E+1/6.0E+1)+t111*t136+t85*t119*t136-t84*t95*t105*t136*(3.0/5.0)-(t89*t96*t108*t136)/5.0+t85*t110*t138*t140+t84*t112*t138*t140+t89*t113*t138*t140
    dqLdq(4) = t100*t136*(-1.0/1.0E+1)-t107*t136*(7.0/6.0E+1)+t89*t127*t136-t85*t90*t101*t136*(6.0/5.0)+t85*t110*t138*(t115-t127)+t84*t112*t138*(t115-t127)+t89*t113*t138*(t115-t127)
    dqLdq(5) = (t107*t136)/3.0E+1-(t89*t93*t108*t136)/5.0+(t89*t93*t109*t138)/5.0E+1+t85*t93*t100*t108*t138*(3.0/2.5E+1)+t84*t93*t104*t108*t138*(3.0/5.0E+1)
    dqRdq(1) = (t104*t135)/3.0E+1-(t88*t92*t105*t135)/5.0+(t88*t92*t106*t137)/5.0E+1+t86*t92*t100*t105*t137*(3.0/2.5E+1)+t87*t92*t105*t107*t137*(3.0/5.0E+1)
    dqRdq(2) = t100*t135*(-1.0/1.0E+1)-t104*t135*(7.0/6.0E+1)+t88*t125*t135-t86*t91*t101*t135*(6.0/5.0)+t86*t110*t137*(t116-t125)+t88*t111*t137*(t116-t125)+t87*t114*t137*(t116-t125)
    dqRdq(3) = (t100*t135)/2.0+t104*t135*(1.1E+1/6.0E+1)+t113*t135+t86*t119*t135-(t88*t95*t105*t135)/5.0-t87*t96*t108*t135*(3.0/5.0)+t86*t110*t137*t139+t88*t111*t137*t139+t87*t114*t137*t139
    dqRdq(4) = (t100*t135)/5.0+(t107*t135)/4.0+t87*t128*t135-t86*t90*t101*t135*(6.0/5.0)+t86*t110*t137*(t115-t128)+t88*t111*t137*(t115-t128)+t87*t114*t137*(t115-t128)
    dqRdq(5) = t107*t135*(-1.0/2.0E+1)-t87*t93*t108*t135*(3.0/5.0)+t87*t93*t109*t137*(9.0/5.0E+1)+t86*t93*t100*t108*t137*(9.0/2.5E+1)+t88*t93*t104*t108*t137*(3.0/5.0E+1)
    
  end subroutine WENO5_JS_jab
  
  subroutine WENO5_M_jab(this,dqLdq,dqRdq,q)
    class(recon_type), intent(in) :: this
    real(8), dimension(this%nc), intent(inout) :: dqLdq,dqRdq
    real(8), dimension(this%nc), intent(in   ) :: q
    
    real(8) :: q1,q2,q3,q4,q5
    
    real(8), parameter :: c1      = 0.1
    real(8), parameter :: c2      = 0.6
    real(8), parameter :: c3      = 0.3
    
    real(8) :: &
    t2    ,&
    t3    ,&
    t4    ,&
    t5    ,&
    t6    ,&
    t7    ,&
    t10   ,&
    t14   ,&
    t15   ,&
    t16   ,&
    t17   ,&
    t18   ,&
    t19   ,&
    t20   ,&
    t21   ,&
    t22   ,&
    t23   ,&
    t24   ,&
    t25   ,&
    t26   ,&
    t27   ,&
    t28   ,&
    t29   ,&
    t34   ,&
    t36   ,&
    t37   ,&
    t39   ,&
    t40   ,&
    t41   ,&
    t43   ,&
    t44   ,&
    t45   ,&
    t46   ,&
    t47   ,&
    t48   ,&
    t49   ,&
    t50   ,&
    t51   ,&
    t52   ,&
    t53   ,&
    t54   ,&
    t55   ,&
    t56   ,&
    t57   ,&
    t58   ,&
    t8   ,&
    t9   ,&
    t11   ,&
    t12   ,&
    t13   ,&
    t30   ,&
    t31   ,&
    t32   ,&
    t33   ,&
    t35   ,&
    t38   ,&
    t42   ,&
    t62   ,&
    t63   ,&
    t64   ,&
    t65   ,&
    t66   ,&
    t67   ,&
    t68   ,&
    t69   ,&
    t59   ,&
    t60   ,&
    t61   ,&
    t70   ,&
    t71   ,&
    t72   ,&
    t84   ,&
    t85   ,&
    t86   ,&
    t87   ,&
    t88   ,&
    t89   ,&
    t90   ,&
    t91   ,&
    t92   ,&
    t93   ,&
    t94   ,&
    t95   ,&
    t96   ,&
    t97   ,&
    t98   ,&
    t73   ,&
    t74   ,&
    t75   ,&
    t76   ,&
    t77   ,&
    t78   ,&
    t79   ,&
    t80   ,&
    t81   ,&
    t82   ,&
    t83   ,&
    t99   ,&
    t104  ,&
    t105  ,&
    t100  ,&
    t101  ,&
    t103  ,&
    t106  ,&
    t107  ,&
    t109  ,&
    t110  ,&
    t111  ,&
    t112  ,&
    t114  ,&
    t115  ,&
    t102  ,&
    t108  ,&
    t113  ,&
    t116  ,&
    t117  ,&
    t118  ,&
    t119  ,&
    t120  ,&
    t121  ,&
    t122  ,&
    t125  ,&
    t127  ,&
    t128  ,&
    t129  ,&
    t130  ,&
    t131  ,&
    t132  ,&
    t133  ,&
    t134  ,&
    t123  ,&
    t124  ,&
    t126  ,&
    t139  ,&
    t140  ,&
    t135  ,&
    t136  ,&
    t137  ,&
    t138  ,&
    t141  ,&
    t142  ,&
    t193  ,&
    t194  ,&
    t143  ,&
    t144  ,&
    t145  ,&
    t146  ,&
    t147  ,&
    t148  ,&
    t151  ,&
    t152  ,&
    t155  ,&
    t156  ,&
    t157  ,&
    t158  ,&
    t159  ,&
    t161  ,&
    t167  ,&
    t168  ,&
    t189  ,&
    t190  ,&
    t191  ,&
    t192  ,&
    t195  ,&
    t196  ,&
    t201  ,&
    t202  ,&
    t203  ,&
    t204  ,&
    t211  ,&
    t212  ,&
    t215  ,&
    t216  ,&
    t217  ,&
    t218  ,&
    t219  ,&
    t220  ,&
    t222  ,&
    t223  ,&
    t224  ,&
    t226  ,&
    t227  ,&
    t228  ,&
    t230  ,&
    t231  ,&
    t232  ,&
    t234  ,&
    t235  ,&
    t236  ,&
    t243  ,&
    t247  ,&
    t250  ,&
    t255  ,&
    t258  ,&
    t260  ,&
    t149  ,&
    t150  ,&
    t153  ,&
    t154  ,&
    t160  ,&
    t162  ,&
    t163  ,&
    t164  ,&
    t165  ,&
    t166  ,&
    t169  ,&
    t170  ,&
    t171  ,&
    t172  ,&
    t173  ,&
    t174  ,&
    t175  ,&
    t176  ,&
    t197  ,&
    t198  ,&
    t199  ,&
    t200  ,&
    t209  ,&
    t210  ,&
    t221  ,&
    t225  ,&
    t229  ,&
    t233  ,&
    t237  ,&
    t238  ,&
    t239  ,&
    t240  ,&
    t241  ,&
    t242  ,&
    t244  ,&
    t245  ,&
    t246  ,&
    t248  ,&
    t249  ,&
    t251  ,&
    t252  ,&
    t253  ,&
    t254  ,&
    t256  ,&
    t257  ,&
    t259  ,&
    t261  ,&
    t262  ,&
    t263  ,&
    t264  ,&
    t266  ,&
    t267  ,&
    t269  ,&
    t270  ,&
    t271  ,&
    t272  ,&
    t273  ,&
    t274  ,&
    t279  ,&
    t280  ,&
    t283  ,&
    t284  ,&
    t285  ,&
    t286  ,&
    t287  ,&
    t288  ,&
    t293  ,&
    t294  ,&
    t295  ,&
    t296  ,&
    t298  ,&
    t300  ,&
    t301  ,&
    t302  ,&
    t303  ,&
    t304  ,&
    t307  ,&
    t308  ,&
    t309  ,&
    t310  ,&
    t311  ,&
    t312  ,&
    t313  ,&
    t314  ,&
    t315  ,&
    t316  ,&
    t317  ,&
    t318  ,&
    t319  ,&
    t320  ,&
    t321  ,&
    t322  ,&
    t323  ,&
    t324  ,&
    t325  ,&
    t326  ,&
    t327  ,&
    t328  ,&
    t329  ,&
    t330  ,&
    t331  ,&
    t332  ,&
    t333  ,&
    t334  ,&
    t335  ,&
    t336  ,&
    t337  ,&
    t338  ,&
    t339  ,&
    t340  ,&
    t341  ,&
    t342  ,&
    t343  ,&
    t344  ,&
    t345  ,&
    t346  ,&
    t347  ,&
    t348  ,&
    t349  ,&
    t350  ,&
    t357  ,&
    t358  ,&
    t359  ,&
    t360  ,&
    t361  ,&
    t362  ,&
    t365  ,&
    t366  ,&
    t367  ,&
    t368  ,&
    t369  ,&
    t370  ,&
    t372  ,&
    t374  ,&
    t375  ,&
    t376  ,&
    t379  ,&
    t380  ,&
    t177  ,&
    t178  ,&
    t181  ,&
    t183  ,&
    t185  ,&
    t187  ,&
    t205  ,&
    t206  ,&
    t207  ,&
    t208  ,&
    t213  ,&
    t214  ,&
    t265  ,&
    t268  ,&
    t275  ,&
    t276  ,&
    t277  ,&
    t278  ,&
    t281  ,&
    t282  ,&
    t289  ,&
    t290  ,&
    t291  ,&
    t292  ,&
    t297  ,&
    t299  ,&
    t305  ,&
    t306  ,&
    t351  ,&
    t352  ,&
    t353  ,&
    t354  ,&
    t355  ,&
    t356  ,&
    t363  ,&
    t364  ,&
    t371  ,&
    t373  ,&
    t377  ,&
    t378  ,&
    t393  ,&
    t394  ,&
    t395  ,&
    t396  ,&
    t397  ,&
    t398  ,&
    t399  ,&
    t400  ,&
    t401  ,&
    t402  ,&
    t403  ,&
    t404  ,&
    t405  ,&
    t406  ,&
    t411  ,&
    t412  ,&
    t413  ,&
    t414  ,&
    t531  ,&
    t532  ,&
    t533  ,&
    t534  ,&
    t179  ,&
    t180  ,&
    t182  ,&
    t184  ,&
    t186  ,&
    t188  ,&
    t381  ,&
    t382  ,&
    t383  ,&
    t384  ,&
    t385  ,&
    t386  ,&
    t387  ,&
    t388  ,&
    t389  ,&
    t390  ,&
    t391  ,&
    t392  ,&
    t407  ,&
    t408  ,&
    t409  ,&
    t410  ,&
    t415  ,&
    t416  ,&
    t417  ,&
    t418  ,&
    t419  ,&
    t420  ,&
    t421  ,&
    t422  ,&
    t423  ,&
    t424  ,&
    t425  ,&
    t426  ,&
    t427  ,&
    t428  ,&
    t429  ,&
    t430  ,&
    t431  ,&
    t432  ,&
    t433  ,&
    t434  ,&
    t436  ,&
    t438  ,&
    t439  ,&
    t440  ,&
    t441  ,&
    t442  ,&
    t443  ,&
    t444  ,&
    t447  ,&
    t448  ,&
    t449  ,&
    t450  ,&
    t452  ,&
    t454  ,&
    t461  ,&
    t462  ,&
    t465  ,&
    t466  ,&
    t473  ,&
    t474  ,&
    t475  ,&
    t476  ,&
    t477  ,&
    t478  ,&
    t479  ,&
    t480  ,&
    t481  ,&
    t482  ,&
    t483  ,&
    t484  ,&
    t485  ,&
    t486  ,&
    t487  ,&
    t488  ,&
    t489  ,&
    t490  ,&
    t491  ,&
    t492  ,&
    t493  ,&
    t494  ,&
    t495  ,&
    t496  ,&
    t499  ,&
    t500  ,&
    t501  ,&
    t502  ,&
    t503  ,&
    t504  ,&
    t505  ,&
    t506  ,&
    t509  ,&
    t510  ,&
    t511  ,&
    t512  ,&
    t513  ,&
    t514  ,&
    t515  ,&
    t516  ,&
    t517  ,&
    t518  ,&
    t523  ,&
    t524  ,&
    t525  ,&
    t526  ,&
    t527  ,&
    t528  ,&
    t529  ,&
    t530  ,&
    t539  ,&
    t540  ,&
    t541  ,&
    t542  ,&
    t543  ,&
    t544  ,&
    t565  ,&
    t566  ,&
    t567  ,&
    t568  ,&
    t569  ,&
    t570  ,&
    t571  ,&
    t572  ,&
    t573  ,&
    t574  ,&
    t575  ,&
    t576  ,&
    t577  ,&
    t578  ,&
    t579  ,&
    t580  ,&
    t435  ,&
    t437  ,&
    t445  ,&
    t446  ,&
    t451  ,&
    t453  ,&
    t455  ,&
    t456  ,&
    t457  ,&
    t458  ,&
    t459  ,&
    t460  ,&
    t463  ,&
    t464  ,&
    t469  ,&
    t470  ,&
    t497  ,&
    t498  ,&
    t507  ,&
    t508  ,&
    t519  ,&
    t520  ,&
    t521  ,&
    t522  ,&
    t535  ,&
    t536  ,&
    t537  ,&
    t538  ,&
    t545  ,&
    t546  ,&
    t547  ,&
    t548  ,&
    t549  ,&
    t550  ,&
    t551  ,&
    t552  ,&
    t553  ,&
    t554  ,&
    t555  ,&
    t556  ,&
    t557  ,&
    t558  ,&
    t559  ,&
    t560  ,&
    t561  ,&
    t562  ,&
    t563  ,&
    t564  ,&
    t581  ,&
    t582  ,&
    t583  ,&
    t584  ,&
    t585  ,&
    t586  ,&
    t467  ,&
    t468  ,&
    t471  ,&
    t472  ,&
    t587  ,&
    t588  ,&
    t595  ,&
    t596  ,&
    t597  ,&
    t598  ,&
    t599  ,&
    t600  ,&
    t589  ,&
    t590  ,&
    t591  ,&
    t592  ,&
    t593  ,&
    t594
    
    q1 = q(1)
    q2 = q(2)
    q3 = q(3)
    q4 = q(4)
    q5 = q(5)
    
    t2 = q2*2.0
    t3 = q3*2.0
    t4 = q2*4.0
    t5 = q3*3.0
    t6 = q4*2.0
    t7 = q4*4.0
    t10 = -q4
    t14 = q1/3.0
    t15 = q2/3.0
    t16 = q3/3.0
    t17 = q1/6.0
    t18 = q4/3.0
    t19 = q2/6.0
    t20 = q5/3.0
    t21 = q2*(5.0/3.0)
    t22 = q4/6.0
    t23 = q4*(5.0/3.0)
    t24 = q5/6.0
    t25 = q2*(5.0/6.0)
    t26 = q3*(5.0/6.0)
    t27 = q2*(7.0/6.0)
    t28 = q4*(5.0/6.0)
    t29 = q4*(7.0/6.0)
    t34 = q1*(8.0/3.0)
    t36 = q2*(8.0/3.0)
    t37 = q1*(1.1E+1/3.0)
    t39 = q4*(8.0/3.0)
    t40 = q5*(8.0/3.0)
    t41 = q3*(1.1E+1/3.0)
    t43 = q2*(1.3E+1/3.0)
    t44 = q3*(1.3E+1/3.0)
    t45 = q5*(1.1E+1/3.0)
    t46 = q3*(1.1E+1/6.0)
    t47 = q4*(1.3E+1/3.0)
    t48 = q1*(1.9E+1/3.0)
    t49 = q2*(1.9E+1/3.0)
    t50 = q3*(2.0E+1/3.0)
    t51 = q4*(1.9E+1/3.0)
    t52 = q5*(1.9E+1/3.0)
    t53 = q3*(2.6E+1/3.0)
    t54 = q2*(3.1E+1/3.0)
    t55 = q3*(3.1E+1/3.0)
    t56 = q4*(3.1E+1/3.0)
    t57 = q2*(5.0E+1/3.0)
    t58 = q4*(5.0E+1/3.0)
    t8 = -t2
    t9 = -t3
    t11 = -t4
    t12 = -t6
    t13 = -t7
    t30 = q2+t10
    t31 = -t17
    t32 = -t19
    t33 = -t22
    t35 = -t24
    t38 = -t27
    t42 = -t29
    t62 = -t44
    t63 = -t49
    t64 = -t51
    t65 = -t53
    t66 = -t54
    t67 = -t56
    t68 = -t57
    t69 = -t58
    t59 = q1+q3+t8
    t60 = q2+q4+t9
    t61 = q3+q5+t12
    t70 = t30*t30
    t71 = q1+t5+t11
    t72 = q5+t5+t13
    t84 = t16+t25+t31
    t85 = t15+t26+t33
    t86 = t18+t26+t32
    t87 = t16+t28+t35
    t88 = t14+t38+t46
    t89 = t20+t42+t46
    t90 = t21+t39+t62
    t91 = t23+t36+t62
    t92 = t34+t41+t63
    t93 = t40+t41+t64
    t94 = t43+t47+t65
    t95 = t37+t50+t66
    t96 = t45+t50+t67
    t97 = t48+t55+t68
    t98 = t52+t55+t69
    t73 = t59*t59
    t74 = t60*t60
    t75 = t61*t61
    t76 = t71*t71
    t77 = t72*t72
    t78 = t70/4.0
    t79 = t76/4.0
    t80 = t77/4.0
    t81 = t73*(1.3E+1/1.2E+1)
    t82 = t74*(1.3E+1/1.2E+1)
    t83 = t75*(1.3E+1/1.2E+1)
    t99 = this%eps+t78+t82
    t104 = this%eps+t79+t81
    t105 = this%eps+t80+t83
    t100 = 1.0/(t99*t99)
    t101 = 1.0/(t99*t99*t99)
    t103 = 1.0/(t99*t99*t99*t99*t99)
    t106 = 1.0/(t104*t104)
    t107 = 1.0/(t104*t104*t104)
    t109 = 1.0/(t104*t104*t104*t104*t104)
    t110 = 1.0/(t104*t104*t104*t104*t104*t104*t104)
    t111 = 1.0/(t105*t105)
    t112 = 1.0/(t105*t105*t105)
    t114 = 1.0/(t105*t105*t105*t105*t105)
    t115 = 1.0/(t105*t105*t105*t105*t105*t105*t105)
    t102 = t100*t100
    t108 = t106*t106
    t113 = t111*t111
    t116 = t100*(3.0/5.0)
    t117 = t106/1.0E+1
    t118 = t106*(3.0/1.0E+1)
    t119 = t111/1.0E+1
    t120 = t111*(3.0/1.0E+1)
    t121 = t90*t101*(6.0/5.0)
    t122 = t91*t101*(6.0/5.0)
    t125 = t94*t101*(6.0/5.0)
    t127 = (t95*t107)/5.0
    t128 = t95*t107*(3.0/5.0)
    t129 = (t96*t112)/5.0
    t130 = t96*t112*(3.0/5.0)
    t131 = (t97*t107)/5.0
    t132 = t97*t107*(3.0/5.0)
    t133 = (t98*t112)/5.0
    t134 = t98*t112*(3.0/5.0)
    t123 = -t121
    t124 = -t122
    t126 = -t125
    t139 = t116+t117+t120
    t140 = t116+t118+t119
    t135 = t124+t131
    t136 = t124+t132
    t137 = t123+t133
    t138 = t123+t134
    t141 = 1.0/t139
    t142 = 1.0/t140
    t193 = t126+t127+t130
    t194 = t126+t128+t129
    t143 = t141*t141
    t144 = t142*t142
    t145 = t141*t141*t141
    t146 = t142*t142*t142
    t147 = t100*t141*(3.0/2.5E+1)
    t148 = t100*t142*(3.0/2.5E+1)
    t151 = t100*t141*(2.7E+1/2.5E+1)
    t152 = t100*t142*(2.7E+1/2.5E+1)
    t155 = t106*t141*(2.0/2.5E+1)
    t156 = t106*t142*(3.0/2.5E+1)
    t157 = t111*t142*(2.0/2.5E+1)
    t158 = t111*t141*(3.0/2.5E+1)
    t159 = t106*t141*(3.0/1.0E+2)
    t161 = t111*t142*(3.0/1.0E+2)
    t167 = t106*t142*(2.7E+1/1.0E+2)
    t168 = t111*t141*(2.7E+1/1.0E+2)
    t189 = t90*t101*t141*(6.0/2.5E+1)
    t190 = t91*t101*t141*(6.0/2.5E+1)
    t191 = t90*t101*t142*(6.0/2.5E+1)
    t192 = t91*t101*t142*(6.0/2.5E+1)
    t195 = t94*t101*t141*(6.0/2.5E+1)
    t196 = t94*t101*t142*(6.0/2.5E+1)
    t201 = t90*t101*t141*(5.4E+1/2.5E+1)
    t202 = t91*t101*t141*(5.4E+1/2.5E+1)
    t203 = t90*t101*t142*(5.4E+1/2.5E+1)
    t204 = t91*t101*t142*(5.4E+1/2.5E+1)
    t211 = t94*t101*t141*(5.4E+1/2.5E+1)
    t212 = t94*t101*t142*(5.4E+1/2.5E+1)
    t215 = t92*t107*t141*(4.0/2.5E+1)
    t216 = t92*t107*t142*(6.0/2.5E+1)
    t217 = t93*t112*t142*(4.0/2.5E+1)
    t218 = t93*t112*t141*(6.0/2.5E+1)
    t219 = t92*t107*t141*(3.0/5.0E+1)
    t220 = t95*t107*t141*(4.0/2.5E+1)
    t222 = t95*t107*t142*(6.0/2.5E+1)
    t223 = t93*t112*t142*(3.0/5.0E+1)
    t224 = t96*t112*t142*(4.0/2.5E+1)
    t226 = t96*t112*t141*(6.0/2.5E+1)
    t227 = t95*t107*t141*(3.0/5.0E+1)
    t228 = t97*t107*t141*(4.0/2.5E+1)
    t230 = t97*t107*t142*(6.0/2.5E+1)
    t231 = t96*t112*t142*(3.0/5.0E+1)
    t232 = t98*t112*t142*(4.0/2.5E+1)
    t234 = t98*t112*t141*(6.0/2.5E+1)
    t235 = t97*t107*t141*(3.0/5.0E+1)
    t236 = t98*t112*t142*(3.0/5.0E+1)
    t243 = t92*t107*t142*(2.7E+1/5.0E+1)
    t247 = t93*t112*t141*(2.7E+1/5.0E+1)
    t250 = t95*t107*t142*(2.7E+1/5.0E+1)
    t255 = t96*t112*t141*(2.7E+1/5.0E+1)
    t258 = t97*t107*t142*(2.7E+1/5.0E+1)
    t260 = t98*t112*t141*(2.7E+1/5.0E+1)
    t149 = t102*t143*(9.0/2.5E+1)
    t150 = t102*t144*(9.0/2.5E+1)
    t153 = -t151
    t154 = -t152
    t160 = (t108*t143)/1.0E+2
    t162 = (t113*t144)/1.0E+2
    t163 = -t159
    t164 = t108*t144*(9.0/1.0E+2)
    t165 = -t161
    t166 = t113*t143*(9.0/1.0E+2)
    t169 = -t167
    t170 = -t168
    t171 = t147-9.0/2.5E+1
    t172 = t148-9.0/2.5E+1
    t173 = t155+1.0/1.0E+2
    t174 = t157+1.0/1.0E+2
    t175 = t156+9.0/1.0E+2
    t176 = t158+9.0/1.0E+2
    t197 = t90*t103*t143*(3.6E+1/2.5E+1)
    t198 = t91*t103*t143*(3.6E+1/2.5E+1)
    t199 = t90*t103*t144*(3.6E+1/2.5E+1)
    t200 = t91*t103*t144*(3.6E+1/2.5E+1)
    t209 = t94*t103*t143*(3.6E+1/2.5E+1)
    t210 = t94*t103*t144*(3.6E+1/2.5E+1)
    t221 = (t95*t109*t143)/2.5E+1
    t225 = (t96*t114*t144)/2.5E+1
    t229 = (t97*t109*t143)/2.5E+1
    t233 = (t98*t114*t144)/2.5E+1
    t237 = -t220
    t238 = -t222
    t239 = t95*t109*t144*(9.0/2.5E+1)
    t240 = -t224
    t241 = -t226
    t242 = -t227
    t244 = t96*t114*t143*(9.0/2.5E+1)
    t245 = -t228
    t246 = -t230
    t248 = -t231
    t249 = t97*t109*t144*(9.0/2.5E+1)
    t251 = -t232
    t252 = -t234
    t253 = -t235
    t254 = t98*t114*t143*(9.0/2.5E+1)
    t256 = t92*t109*t143*(2.0/1.25E+2)
    t257 = -t236
    t259 = t93*t114*t144*(2.0/1.25E+2)
    t261 = (t92*t110*t145)/2.5E+2
    t262 = (t93*t115*t146)/2.5E+2
    t263 = -t250
    t264 = -t255
    t266 = t92*t109*t144*(9.0/1.25E+2)
    t267 = -t258
    t269 = t93*t114*t143*(9.0/1.25E+2)
    t270 = -t260
    t271 = t92*t110*t146*(2.7E+1/2.5E+2)
    t272 = t93*t115*t145*(2.7E+1/2.5E+2)
    t273 = t92*t109*t143*(2.3E+1/5.0E+2)
    t274 = t93*t114*t144*(2.3E+1/5.0E+2)
    t279 = t92*t109*t144*(2.61E+2/5.0E+2)
    t280 = t93*t114*t143*(2.61E+2/5.0E+2)
    t283 = t92*t102*t107*t145*(1.8E+1/1.25E+2)
    t284 = t92*t100*t107*t143*(2.7E+1/1.25E+2)
    t285 = t93*t102*t112*t146*(1.8E+1/1.25E+2)
    t286 = t93*t100*t112*t144*(2.7E+1/1.25E+2)
    t287 = t92*t102*t107*t146*(5.4E+1/1.25E+2)
    t288 = t93*t102*t112*t145*(5.4E+1/1.25E+2)
    t293 = t92*t100*t107*t144*(8.1E+1/1.25E+2)
    t294 = t93*t100*t112*t143*(8.1E+1/1.25E+2)
    t295 = t92*t107*t113*t146*(3.0/2.5E+2)
    t296 = t93*t108*t112*t145*(3.0/2.5E+2)
    t298 = t92*t107*t113*t145*(9.0/2.5E+2)
    t300 = t93*t108*t112*t146*(9.0/2.5E+2)
    t301 = t92*t107*t111*t144*(9.0/5.0E+2)
    t302 = t93*t106*t112*t143*(9.0/5.0E+2)
    t303 = t92*t107*t111*t143*(2.7E+1/5.0E+2)
    t304 = t93*t106*t112*t144*(2.7E+1/5.0E+2)
    t307 = t100*t143*(t122-t131)*(-3.0/2.5E+1)
    t308 = t100*t144*(t122-t132)*(-3.0/2.5E+1)
    t309 = t100*t144*(t121-t133)*(-3.0/2.5E+1)
    t310 = t100*t143*(t121-t134)*(-3.0/2.5E+1)
    t311 = t102*t145*(t122-t131)*(-1.8E+1/2.5E+1)
    t312 = t102*t146*(t122-t132)*(-1.8E+1/2.5E+1)
    t313 = t100*t143*(t122-t131)*(-2.7E+1/2.5E+1)
    t314 = t100*t144*(t122-t132)*(-2.7E+1/2.5E+1)
    t315 = t102*t146*(t121-t133)*(-1.8E+1/2.5E+1)
    t316 = t102*t145*(t121-t134)*(-1.8E+1/2.5E+1)
    t317 = t100*t144*(t121-t133)*(-2.7E+1/2.5E+1)
    t318 = t100*t143*(t121-t134)*(-2.7E+1/2.5E+1)
    t319 = t102*t145*(t122-t131)*(1.8E+1/2.5E+1)
    t320 = t102*t146*(t122-t132)*(1.8E+1/2.5E+1)
    t321 = t102*t146*(t121-t133)*(1.8E+1/2.5E+1)
    t322 = t102*t145*(t121-t134)*(1.8E+1/2.5E+1)
    t323 = t106*t143*(t122-t131)*(-2.0/2.5E+1)
    t324 = t106*t144*(t122-t132)*(-3.0/2.5E+1)
    t325 = t108*t145*(t122-t131)*(-1.0/5.0E+1)
    t326 = t111*t144*(t121-t133)*(-2.0/2.5E+1)
    t327 = t111*t143*(t121-t134)*(-3.0/2.5E+1)
    t328 = t113*t146*(t122-t132)*(-1.0/5.0E+1)
    t329 = t108*t145*(t121-t134)*(-1.0/5.0E+1)
    t330 = t113*t146*(t121-t133)*(-1.0/5.0E+1)
    t331 = (t108*t145*(t122-t131))/5.0E+1
    t332 = t108*t146*(t122-t132)*(-9.0/5.0E+1)
    t333 = (t113*t146*(t122-t132))/5.0E+1
    t334 = t113*t145*(t122-t131)*(-9.0/5.0E+1)
    t335 = (t108*t145*(t121-t134))/5.0E+1
    t336 = t108*t146*(t121-t133)*(-9.0/5.0E+1)
    t337 = (t113*t146*(t121-t133))/5.0E+1
    t338 = t113*t145*(t121-t134)*(-9.0/5.0E+1)
    t339 = t106*t143*(t122-t131)*(-3.0/1.0E+2)
    t340 = t111*t144*(t122-t132)*(-3.0/1.0E+2)
    t341 = t106*t143*(t121-t134)*(-3.0/1.0E+2)
    t342 = t111*t144*(t121-t133)*(-3.0/1.0E+2)
    t343 = t108*t146*(t122-t132)*(9.0/5.0E+1)
    t344 = t113*t145*(t122-t131)*(9.0/5.0E+1)
    t345 = t108*t146*(t121-t133)*(9.0/5.0E+1)
    t346 = t113*t145*(t121-t134)*(9.0/5.0E+1)
    t347 = t106*t144*(t122-t132)*(-2.7E+1/1.0E+2)
    t348 = t111*t143*(t122-t131)*(-2.7E+1/1.0E+2)
    t349 = t106*t144*(t121-t133)*(-2.7E+1/1.0E+2)
    t350 = t111*t143*(t121-t134)*(-2.7E+1/1.0E+2)
    t357 = t100*t143*t193*(3.0/2.5E+1)
    t358 = t100*t144*t194*(3.0/2.5E+1)
    t359 = t102*t145*t193*(1.8E+1/2.5E+1)
    t360 = t102*t146*t194*(1.8E+1/2.5E+1)
    t361 = t100*t143*t193*(2.7E+1/2.5E+1)
    t362 = t100*t144*t194*(2.7E+1/2.5E+1)
    t365 = t106*t143*t193*(2.0/2.5E+1)
    t366 = t106*t144*t194*(3.0/2.5E+1)
    t367 = t111*t144*t194*(2.0/2.5E+1)
    t368 = t111*t143*t193*(3.0/2.5E+1)
    t369 = (t108*t145*t193)/5.0E+1
    t370 = (t113*t146*t194)/5.0E+1
    t372 = t108*t146*t194*(9.0/5.0E+1)
    t374 = t113*t145*t193*(9.0/5.0E+1)
    t375 = t106*t143*t193*(3.0/1.0E+2)
    t376 = t111*t144*t194*(3.0/1.0E+2)
    t379 = t106*t144*t194*(2.7E+1/1.0E+2)
    t380 = t111*t143*t193*(2.7E+1/1.0E+2)
    t177 = 1.0/t171
    t178 = 1.0/t172
    t181 = 1.0/t173
    t183 = 1.0/t174
    t185 = 1.0/t175
    t187 = 1.0/t176
    t205 = -t197
    t206 = -t198
    t207 = -t199
    t208 = -t200
    t213 = -t209
    t214 = -t210
    t265 = -t256
    t268 = -t259
    t275 = -t266
    t276 = -t269
    t277 = -t273
    t278 = -t274
    t281 = -t279
    t282 = -t280
    t289 = -t283
    t290 = -t285
    t291 = -t287
    t292 = -t288
    t297 = -t295
    t299 = -t296
    t305 = -t298
    t306 = -t300
    t351 = t149+t153+2.4E+1/2.5E+1
    t352 = t150+t154+2.4E+1/2.5E+1
    t353 = t160+t163+1.1E+1/1.0E+2
    t354 = t162+t165+1.1E+1/1.0E+2
    t355 = t164+t169+3.9E+1/1.0E+2
    t356 = t166+t170+3.9E+1/1.0E+2
    t363 = -t359
    t364 = -t360
    t371 = -t369
    t373 = -t370
    t377 = -t372
    t378 = -t374
    t393 = t190+t307
    t394 = t192+t308
    t395 = t191+t309
    t396 = t189+t310
    t397 = t245+t323
    t398 = t246+t324
    t399 = t251+t326
    t400 = t252+t327
    t401 = t195+t357
    t402 = t196+t358
    t403 = t237+t365
    t404 = t238+t366
    t405 = t240+t367
    t406 = t241+t368
    t411 = t333+t340
    t412 = t335+t341
    t413 = t344+t348
    t414 = t345+t349
    t531 = t229+t253+t331+t339
    t532 = t233+t257+t337+t342
    t533 = t249+t267+t343+t347
    t534 = t254+t270+t346+t350
    t179 = t177*t177
    t180 = t178*t178
    t182 = t181*t181
    t184 = t183*t183
    t186 = t185*t185
    t188 = t187*t187
    t381 = t215+t265
    t382 = t217+t268
    t383 = t216+t275
    t384 = t218+t276
    t385 = t284+t289
    t386 = t286+t290
    t387 = t291+t293
    t388 = t292+t294
    t389 = t297+t301
    t390 = t299+t302
    t391 = t303+t305
    t392 = t304+t306
    t407 = t219+t261+t277
    t408 = t223+t262+t278
    t409 = t243+t271+t281
    t410 = t247+t272+t282
    t415 = t116*t141*t177*t351
    t416 = t116*t142*t178*t352
    t417 = t100*t141*t177*t351*(-3.0/5.0)
    t418 = t100*t142*t178*t352*(-3.0/5.0)
    t419 = t117*t141*t181*t353
    t420 = t119*t142*t183*t354
    t421 = t118*t142*t185*t355
    t422 = t120*t141*t187*t356
    t423 = t121*t141*t177*t351
    t424 = t122*t141*t177*t351
    t425 = t121*t142*t178*t352
    t426 = t122*t142*t178*t352
    t427 = t125*t141*t177*t351
    t428 = t125*t142*t178*t352
    t429 = (t92*t107*t141*t181*t353)/5.0
    t430 = t127*t141*t181*t353
    t431 = (t93*t112*t142*t183*t354)/5.0
    t432 = t131*t141*t181*t353
    t433 = t129*t142*t183*t354
    t434 = t133*t142*t183*t354
    t436 = (t92*t109*t143*t181*t353)/5.0E+1
    t438 = (t93*t114*t144*t183*t354)/5.0E+1
    t439 = t92*t107*t142*t185*t355*(3.0/5.0)
    t440 = t128*t142*t185*t355
    t441 = t93*t112*t141*t187*t356*(3.0/5.0)
    t442 = t132*t142*t185*t355
    t443 = t130*t141*t187*t356
    t444 = t134*t141*t187*t356
    t447 = t92*t109*t144*t185*t355*(9.0/5.0E+1)
    t448 = t93*t114*t143*t187*t356*(9.0/5.0E+1)
    t449 = t92*t100*t107*t143*t177*t351*(3.0/2.5E+1)
    t450 = t93*t100*t112*t144*t178*t352*(3.0/2.5E+1)
    t452 = t92*t100*t107*t144*t178*t352*(9.0/2.5E+1)
    t454 = t93*t100*t112*t143*t177*t351*(9.0/2.5E+1)
    t461 = t93*t106*t112*t143*t181*t353*(3.0/5.0E+1)
    t462 = t92*t107*t111*t144*t183*t354*(3.0/5.0E+1)
    t465 = t93*t106*t112*t144*t185*t355*(3.0/5.0E+1)
    t466 = t92*t107*t111*t143*t187*t356*(3.0/5.0E+1)
    t473 = -t116*t141*t177*(t283-t284)
    t474 = -t116*t142*t178*(t285-t286)
    t475 = -t116*t142*t178*(t287-t293)
    t476 = -t116*t141*t177*(t288-t294)
    t477 = t100*t143*t177*t351*(t122-t131)*(-3.0/5.0)
    t478 = t100*t144*t178*t352*(t122-t132)*(-3.0/5.0)
    t479 = t100*t144*t178*t352*(t121-t133)*(-3.0/5.0)
    t480 = t100*t143*t177*t351*(t121-t134)*(-3.0/5.0)
    t481 = -t117*t141*t181*(t296-t302)
    t482 = -t119*t142*t183*(t295-t301)
    t483 = (t106*t141*t181*(t296-t302))/1.0E+1
    t484 = (t111*t142*t183*(t295-t301))/1.0E+1
    t485 = -t118*t142*t185*(t300-t304)
    t486 = -t120*t141*t187*(t298-t303)
    t487 = t106*t142*t185*(t300-t304)*(3.0/1.0E+1)
    t488 = t111*t141*t187*(t298-t303)*(3.0/1.0E+1)
    t489 = t106*t143*t181*t353*(t122-t131)*(-1.0/1.0E+1)
    t490 = t106*t143*t181*t353*(t121-t134)*(-1.0/1.0E+1)
    t491 = t111*t144*t183*t354*(t122-t132)*(-1.0/1.0E+1)
    t492 = t111*t144*t183*t354*(t121-t133)*(-1.0/1.0E+1)
    t493 = t117*t143*t181*t353*(t122-t131)
    t494 = t117*t143*t181*t353*(t121-t134)
    t495 = t119*t144*t183*t354*(t122-t132)
    t496 = t119*t144*t183*t354*(t121-t133)
    t499 = t106*t144*t185*t355*(t122-t132)*(-3.0/1.0E+1)
    t500 = t106*t144*t185*t355*(t121-t133)*(-3.0/1.0E+1)
    t501 = t111*t143*t187*t356*(t122-t131)*(-3.0/1.0E+1)
    t502 = t111*t143*t187*t356*(t121-t134)*(-3.0/1.0E+1)
    t503 = t118*t144*t185*t355*(t122-t132)
    t504 = t118*t144*t185*t355*(t121-t133)
    t505 = t120*t143*t187*t356*(t122-t131)
    t506 = t120*t143*t187*t356*(t121-t134)
    t509 = t116*t143*t177*t193*t351
    t510 = t116*t144*t178*t194*t352
    t511 = t117*t143*t181*t193*t353
    t512 = t119*t144*t183*t194*t354
    t513 = t106*t143*t181*t193*t353*(-1.0/1.0E+1)
    t514 = t111*t144*t183*t194*t354*(-1.0/1.0E+1)
    t515 = t118*t144*t185*t194*t355
    t516 = t120*t143*t187*t193*t356
    t517 = t106*t144*t185*t194*t355*(-3.0/1.0E+1)
    t518 = t111*t143*t187*t193*t356*(-3.0/1.0E+1)
    t523 = t106*t141*t181*(t329+t106*t143*(t121-t134)*(3.0/1.0E+2))*(-1.0/1.0E+1)
    t524 = t111*t142*t183*(t328+t111*t144*(t122-t132)*(3.0/1.0E+2))*(-1.0/1.0E+1)
    t525 = t106*t142*t185*(t336+t106*t144*(t121-t133)*(2.7E+1/1.0E+2))*(-3.0/1.0E+1)
    t526 = t111*t141*t187*(t334+t111*t143*(t122-t131)*(2.7E+1/1.0E+2))*(-3.0/1.0E+1)
    t527 = t202+t206+t313+t319
    t528 = t204+t208+t314+t320
    t529 = t203+t207+t317+t321
    t530 = t201+t205+t318+t322
    t539 = t211+t213+t361+t363
    t540 = t212+t214+t362+t364
    t541 = t221+t242+t371+t375
    t542 = t225+t248+t373+t376
    t543 = t239+t263+t377+t379
    t544 = t244+t264+t378+t380
    t565 = -t116*t141*t177*(t198-t202+t311+t100*t143*(t122-t131)*(2.7E+1/2.5E+1))
    t566 = -t116*t142*t178*(t200-t204+t312+t100*t144*(t122-t132)*(2.7E+1/2.5E+1))
    t567 = -t116*t142*t178*(t199-t203+t315+t100*t144*(t121-t133)*(2.7E+1/2.5E+1))
    t568 = -t116*t141*t177*(t197-t201+t316+t100*t143*(t121-t134)*(2.7E+1/2.5E+1))
    t569 = t100*t141*t177*(t198-t202+t311+t100*t143*(t122-t131)*(2.7E+1/2.5E+1))*(3.0/5.0)
    t570 = t100*t142*t178*(t200-t204+t312+t100*t144*(t122-t132)*(2.7E+1/2.5E+1))*(3.0/5.0)
    t571 = t100*t142*t178*(t199-t203+t315+t100*t144*(t121-t133)*(2.7E+1/2.5E+1))*(3.0/5.0)
    t572 = t100*t141*t177*(t197-t201+t316+t100*t143*(t121-t134)*(2.7E+1/2.5E+1))*(3.0/5.0)
    t573 = (t106*t141*t181*t531)/1.0E+1
    t574 = (t111*t142*t183*t532)/1.0E+1
    t575 = t106*t142*t185*t533*(3.0/1.0E+1)
    t576 = t111*t141*t187*t534*(3.0/1.0E+1)
    t577 = -t116*t141*t177*(t209-t211+t359-t361)
    t578 = -t116*t142*t178*(t210-t212+t360-t362)
    t579 = t100*t141*t177*(t209-t211+t359-t361)*(3.0/5.0)
    t580 = t100*t142*t178*(t210-t212+t360-t362)*(3.0/5.0)
    t435 = -t429
    t437 = -t431
    t445 = -t439
    t446 = -t441
    t451 = -t449
    t453 = -t450
    t455 = -t452
    t456 = -t454
    t457 = t92*t102*t107*t145*t179*t351*(9.0/6.25E+2)
    t458 = t92*t102*t107*t146*t180*t352*(2.7E+1/6.25E+2)
    t459 = t93*t102*t112*t146*t180*t352*(9.0/6.25E+2)
    t460 = t93*t102*t112*t145*t179*t351*(2.7E+1/6.25E+2)
    t463 = t93*t108*t112*t145*t182*t353*(3.0/6.25E+2)
    t464 = t92*t107*t113*t146*t184*t354*(3.0/6.25E+2)
    t469 = t93*t108*t112*t146*t186*t355*7.2E-3
    t470 = t92*t107*t113*t145*t188*t356*7.2E-3
    t497 = t108*t145*t182*t353*(t121-t134)*(-1.0/1.25E+2)
    t498 = t113*t146*t184*t354*(t122-t132)*(-1.0/1.25E+2)
    t507 = t108*t146*t186*t355*(t121-t133)*(-9.0/2.5E+2)
    t508 = t113*t145*t188*t356*(t122-t131)*(-9.0/2.5E+2)
    t519 = t117*t141*t181*t407
    t520 = t119*t142*t183*t408
    t521 = t118*t142*t185*t409
    t522 = t120*t141*t187*t410
    t535 = t117*t141*t182*t353*t381
    t536 = t119*t142*t184*t354*t382
    t537 = t118*t142*t186*t355*t383
    t538 = t120*t141*t188*t356*t384
    t545 = t116*t141*t179*t351*t393
    t546 = t116*t142*t180*t352*t394
    t547 = t116*t142*t180*t352*t395
    t548 = t116*t141*t179*t351*t396
    t549 = t100*t141*t179*t351*t393*(-3.0/5.0)
    t550 = t100*t142*t180*t352*t394*(-3.0/5.0)
    t551 = t100*t142*t180*t352*t395*(-3.0/5.0)
    t552 = t100*t141*t179*t351*t396*(-3.0/5.0)
    t553 = t106*t141*t182*t353*(t228+t106*t143*(t122-t131)*(2.0/2.5E+1))*(-1.0/1.0E+1)
    t554 = t111*t142*t184*t354*(t232+t111*t144*(t121-t133)*(2.0/2.5E+1))*(-1.0/1.0E+1)
    t555 = t106*t142*t186*t355*(t230+t106*t144*(t122-t132)*(3.0/2.5E+1))*(-3.0/1.0E+1)
    t556 = t111*t141*t188*t356*(t234+t111*t143*(t121-t134)*(3.0/2.5E+1))*(-3.0/1.0E+1)
    t557 = t116*t141*t179*t351*t401
    t558 = t116*t142*t180*t352*t402
    t559 = t100*t141*t179*t351*t401*(-3.0/5.0)
    t560 = t100*t142*t180*t352*t402*(-3.0/5.0)
    t561 = t106*t141*t182*t353*(t220-t365)*(-1.0/1.0E+1)
    t562 = -t119*t142*t184*t354*(t224-t367)
    t563 = t106*t142*t186*t355*(t222-t366)*(-3.0/1.0E+1)
    t564 = -t120*t141*t188*t356*(t226-t368)
    t581 = (t106*t141*t181*t541)/1.0E+1
    t582 = t119*t142*t183*t542
    t583 = t106*t142*t185*t543*(3.0/1.0E+1)
    t584 = t120*t141*t187*t544
    t585 = t417+t419+t422
    t586 = t418+t420+t421
    t467 = -t463
    t468 = -t464
    t471 = -t469
    t472 = -t470
    t587 = 1.0/t585
    t588 = 1.0/t586
    t595 = t424+t432+t477+t493+t505+t508+t526+t549+t553+t569+t573
    t596 = t425+t434+t479+t496+t504+t507+t525+t551+t554+t571+t574
    t597 = t426+t442+t478+t495+t498+t503+t524+t550+t555+t570+t575
    t598 = t423+t444+t480+t494+t497+t506+t523+t552+t556+t572+t576
    t599 = t427+t430+t443+t509+t513+t518+t559+t561+t564+t579+t581+t584
    t600 = t428+t433+t440+t510+t514+t517+t560+t562+t563+t580+t582+t583
    t589 = t587*t587
    t590 = t588*t588
    t591 = t435+t436+t451+t457+t466+t472+t473+t488+t519+t535
    t592 = t437+t438+t453+t459+t465+t471+t474+t487+t520+t536
    t593 = t445+t447+t455+t458+t462+t468+t475+t484+t521+t537
    t594 = t446+t448+t456+t460+t461+t467+t476+t483+t522+t538
    
    dqLdq(1) = t84*t447*t588+t85*t458*t588+t89*t462*t588+t85*t475*t588+t89*t484*t588+t84*t521*t588+t84*t537*t588+t85*t416*t590*t593-(t106*t142*t185*t355*t588)/2.0E+1-t84*t92*t107*t142*t185*t355*t588*(3.0/5.0)-t84*t106*t142*t185*t355*t590*t593*(3.0/1.0E+1)-(t89*t111*t142*t183*t354*t590*t593)/1.0E+1-t85*t92*t100*t107*t144*t178*t352*t588*(9.0/2.5E+1)-t89*t92*t107*t113*t146*t184*t354*t588*(3.0/6.25E+2)
    dqLdq(2) = t85*t426*t588+t84*t442*t588+t85*t478*t588+t89*t495*t588+t84*t503*t588+t89*t498*t588+t89*t524*t588+t85*t550*t588+t84*t555*t588+t85*t570*t588+t84*t575*t588-t84*t421*t590*(t442+t495+t498+t524-t546+t118*t142*t185*t533+t116*t142*t178*(t200-t204+t312+t100*t144*(t122-t132)*(2.7E+1/2.5E+1))+t91*t101*t142*t178*t352*(6.0/5.0)+t106*t144*t185*t355*(t122-t132)*(3.0/1.0E+1)-t116*t144*t178*t352*(t122-t132)-t118*t142*t186*t355*(t230+t106*t144*(t122-t132)*(3.0/2.5E+1)))-(t100*t142*t178*t352*t588)/5.0+(t106*t142*t185*t355*t588)/4.0+t85*t100*t142*t178*t352*t590*(t442+t495+t498+t524-t546+t118*t142*t185*t533+t116*t142*t178*(t200-t204+t312+t100*t144*(t122-t132)*(2.7E+1/2.5E+1))+t91*t101*t142*t178*t352*(6.0/5.0)+t106*t144*t185*t355*(t122-t132)*(3.0/1.0E+1)-t116*t144*t178*t352*(t122-t132)-t118*t142*t186*t355*(t230+t106*t144*(t122-t132)*(3.0/2.5E+1)))*(3.0/5.0)-(t89*t111*t142*t183*t354*t590*(t442+t495+t498+t524-t546+t118*t142*t185*t533+t116*t142*t178*(t200-t204+t312+t100*t144*(t122-t132)*(2.7E+1/2.5E+1))+t91*t101*t142*t178*t352*(6.0/5.0)+t106*t144*t185*t355*(t122-t132)*(3.0/1.0E+1)-t116*t144*t178*t352*(t122-t132)-t118*t142*t186*t355*(t230+t106*t144*(t122-t132)*(3.0/2.5E+1))))/1.0E+1
    dqLdq(3) = t84*t515*t588+t89*t512*t588+t85*t558*t588+t85*t578*t588+t89*t420*t590*(t440+t510+t514-t515+t560+t562+t580+t582-t126*t142*t178*t352+t118*t142*t185*t543+(t96*t112*t142*t183*t354)/5.0-t118*t142*t186*t355*(t222-t366))-(t100*t142*t178*t352*t588)/2.0+t111*t142*t183*t354*t588*(1.1E+1/6.0E+1)+t117*t142*t185*t355*t588-(t89*t111*t142*t183*t542*t588)/1.0E+1-t84*t118*t142*t185*t543*t588-t85*t100*t142*t178*t352*t590*(t440+t510+t514-t515+t560+t562+t580+t582-t126*t142*t178*t352+t118*t142*t185*t543+(t96*t112*t142*t183*t354)/5.0-t118*t142*t186*t355*(t222-t366))*(3.0/5.0)+t84*t106*t142*t185*t355*t590*(t440+t510+t514-t515+t560+t562+t580+t582-t126*t142*t178*t352+t118*t142*t185*t543+(t96*t112*t142*t183*t354)/5.0-t118*t142*t186*t355*(t222-t366))*(3.0/1.0E+1)-t85*t94*t101*t142*t178*t352*t588*(6.0/5.0)-t84*t95*t107*t142*t185*t355*t588*(3.0/5.0)-(t89*t96*t112*t142*t183*t354*t588)/5.0-t85*t100*t144*t178*t194*t352*t588*(3.0/5.0)+(t89*t111*t142*t184*t354*t588*(t224-t367))/1.0E+1+t84*t118*t142*t186*t355*t588*(t222-t366)
    dqLdq(4) = t85*t425*t588+t89*t434*t588+t85*t479*t588+t89*t496*t588+t84*t504*t588+t84*t507*t588+t84*t525*t588+t85*t551*t588+t89*t554*t588+t85*t571*t588+t89*t574*t588-t84*t421*t590*(t434+t496+t507-t547+t554+t574-t118*t142*t185*(t336+t106*t144*(t121-t133)*(2.7E+1/1.0E+2))+t116*t142*t178*(t199-t203+t315+t100*t144*(t121-t133)*(2.7E+1/2.5E+1))+t90*t101*t142*t178*t352*(6.0/5.0)+t106*t144*t185*t355*(t121-t133)*(3.0/1.0E+1)-t116*t144*t178*t352*(t121-t133))+(t100*t142*t178*t352*t588)/1.0E+1-t111*t142*t183*t354*t588*(7.0/6.0E+1)+t85*t100*t142*t178*t352*t590*(t434+t496+t507-t547+t554+t574-t118*t142*t185*(t336+t106*t144*(t121-t133)*(2.7E+1/1.0E+2))+t116*t142*t178*(t199-t203+t315+t100*t144*(t121-t133)*(2.7E+1/2.5E+1))+t90*t101*t142*t178*t352*(6.0/5.0)+t106*t144*t185*t355*(t121-t133)*(3.0/1.0E+1)-t116*t144*t178*t352*(t121-t133))*(3.0/5.0)-(t89*t111*t142*t183*t354*t590*(t434+t496+t507-t547+t554+t574-t118*t142*t185*(t336+t106*t144*(t121-t133)*(2.7E+1/1.0E+2))+t116*t142*t178*(t199-t203+t315+t100*t144*(t121-t133)*(2.7E+1/2.5E+1))+t90*t101*t142*t178*t352*(6.0/5.0)+t106*t144*t185*t355*(t121-t133)*(3.0/1.0E+1)-t116*t144*t178*t352*(t121-t133)))/1.0E+1
    dqLdq(5) = t89*t438*t588+t85*t459*t588+t84*t465*t588+t85*t474*t588+t84*t487*t588+t89*t520*t588+t89*t536*t588+t85*t416*t590*t592+(t111*t142*t183*t354*t588)/3.0E+1-(t89*t93*t112*t142*t183*t354*t588)/5.0-t84*t106*t142*t185*t355*t590*t592*(3.0/1.0E+1)-(t89*t111*t142*t183*t354*t590*t592)/1.0E+1-t85*t93*t100*t112*t144*t178*t352*t588*(3.0/2.5E+1)-t84*t93*t108*t112*t146*t186*t355*t588*7.2E-3
    dqRdq(1) = t88*t436*t587+t86*t457*t587+t87*t466*t587+t86*t473*t587+t87*t488*t587+t88*t519*t587+t88*t535*t587+t86*t415*t589*t591+(t106*t141*t181*t353*t587)/3.0E+1-(t88*t92*t107*t141*t181*t353*t587)/5.0-(t88*t106*t141*t181*t353*t589*t591)/1.0E+1-t87*t111*t141*t187*t356*t589*t591*(3.0/1.0E+1)-t86*t92*t100*t107*t143*t177*t351*t587*(3.0/2.5E+1)-t87*t92*t107*t113*t145*t188*t356*t587*7.2E-3
    dqRdq(2) = t86*t424*t587+t88*t432*t587+t86*t477*t587+t88*t493*t587+t87*t505*t587+t87*t508*t587+t87*t526*t587+t86*t549*t587+t88*t553*t587+t86*t569*t587+t88*t573*t587-t88*t419*t589*(t432+t505+t508+t526-t545+t117*t141*t181*t531+t116*t141*t177*(t198-t202+t311+t100*t143*(t122-t131)*(2.7E+1/2.5E+1))+t91*t101*t141*t177*t351*(6.0/5.0)+(t106*t143*t181*t353*(t122-t131))/1.0E+1-t116*t143*t177*t351*(t122-t131)-t117*t141*t182*t353*(t228+t106*t143*(t122-t131)*(2.0/2.5E+1)))+(t100*t141*t177*t351*t587)/1.0E+1-t106*t141*t181*t353*t587*(7.0/6.0E+1)+t86*t100*t141*t177*t351*t589*(t432+t505+t508+t526-t545+t117*t141*t181*t531+t116*t141*t177*(t198-t202+t311+t100*t143*(t122-t131)*(2.7E+1/2.5E+1))+t91*t101*t141*t177*t351*(6.0/5.0)+(t106*t143*t181*t353*(t122-t131))/1.0E+1-t116*t143*t177*t351*(t122-t131)-t117*t141*t182*t353*(t228+t106*t143*(t122-t131)*(2.0/2.5E+1)))*(3.0/5.0)-t87*t111*t141*t187*t356*t589*(t432+t505+t508+t526-t545+t117*t141*t181*t531+t116*t141*t177*(t198-t202+t311+t100*t143*(t122-t131)*(2.7E+1/2.5E+1))+t91*t101*t141*t177*t351*(6.0/5.0)+(t106*t143*t181*t353*(t122-t131))/1.0E+1-t116*t143*t177*t351*(t122-t131)-t117*t141*t182*t353*(t228+t106*t143*(t122-t131)*(2.0/2.5E+1)))*(3.0/1.0E+1)
    dqRdq(3) = t88*t511*t587+t87*t516*t587+t86*t557*t587+t86*t577*t587+t87*t422*t589*(t430+t509-t511+t518+t559+t564+t579+t584-t126*t141*t177*t351+t117*t141*t181*t541+t96*t112*t141*t187*t356*(3.0/5.0)-t117*t141*t182*t353*(t220-t365))-(t100*t141*t177*t351*t587)/2.0+t106*t141*t181*t353*t587*(1.1E+1/6.0E+1)+t119*t141*t187*t356*t587-t88*t117*t141*t181*t541*t587-t87*t111*t141*t187*t544*t587*(3.0/1.0E+1)-t86*t100*t141*t177*t351*t589*(t430+t509-t511+t518+t559+t564+t579+t584-t126*t141*t177*t351+t117*t141*t181*t541+t96*t112*t141*t187*t356*(3.0/5.0)-t117*t141*t182*t353*(t220-t365))*(3.0/5.0)+(t88*t106*t141*t181*t353*t589*(t430+t509-t511+t518+t559+t564+t579+t584-t126*t141*t177*t351+t117*t141*t181*t541+t96*t112*t141*t187*t356*(3.0/5.0)-t117*t141*t182*t353*(t220-t365)))/1.0E+1-t86*t94*t101*t141*t177*t351*t587*(6.0/5.0)-(t88*t95*t107*t141*t181*t353*t587)/5.0-t87*t96*t112*t141*t187*t356*t587*(3.0/5.0)-t86*t100*t143*t177*t193*t351*t587*(3.0/5.0)+t88*t117*t141*t182*t353*t587*(t220-t365)+t87*t111*t141*t188*t356*t587*(t226-t368)*(3.0/1.0E+1)
    dqRdq(4) = t86*t423*t587+t87*t444*t587+t86*t480*t587+t88*t494*t587+t88*t497*t587+t87*t506*t587+t88*t523*t587+t86*t552*t587+t87*t556*t587+t86*t572*t587+t87*t576*t587-t88*t419*t589*(t444+t497+t506-t548+t556+t576-t117*t141*t181*(t329+t106*t143*(t121-t134)*(3.0/1.0E+2))+t116*t141*t177*(t197-t201+t316+t100*t143*(t121-t134)*(2.7E+1/2.5E+1))+t90*t101*t141*t177*t351*(6.0/5.0)+(t106*t143*t181*t353*(t121-t134))/1.0E+1-t116*t143*t177*t351*(t121-t134))-(t100*t141*t177*t351*t587)/5.0+(t111*t141*t187*t356*t587)/4.0+t86*t100*t141*t177*t351*t589*(t444+t497+t506-t548+t556+t576-t117*t141*t181*(t329+t106*t143*(t121-t134)*(3.0/1.0E+2))+t116*t141*t177*(t197-t201+t316+t100*t143*(t121-t134)*(2.7E+1/2.5E+1))+t90*t101*t141*t177*t351*(6.0/5.0)+(t106*t143*t181*t353*(t121-t134))/1.0E+1-t116*t143*t177*t351*(t121-t134))*(3.0/5.0)-t87*t111*t141*t187*t356*t589*(t444+t497+t506-t548+t556+t576-t117*t141*t181*(t329+t106*t143*(t121-t134)*(3.0/1.0E+2))+t116*t141*t177*(t197-t201+t316+t100*t143*(t121-t134)*(2.7E+1/2.5E+1))+t90*t101*t141*t177*t351*(6.0/5.0)+(t106*t143*t181*t353*(t121-t134))/1.0E+1-t116*t143*t177*t351*(t121-t134))*(3.0/1.0E+1)
    dqRdq(5) = t87*t448*t587+t86*t460*t587+t88*t461*t587+t86*t476*t587+t88*t483*t587+t87*t522*t587+t87*t538*t587+t86*t415*t589*t594-(t111*t141*t187*t356*t587)/2.0E+1-t87*t93*t112*t141*t187*t356*t587*(3.0/5.0)-(t88*t106*t141*t181*t353*t589*t594)/1.0E+1-t87*t111*t141*t187*t356*t589*t594*(3.0/1.0E+1)-t86*t93*t100*t112*t143*t177*t351*t587*(9.0/2.5E+1)-t88*t93*t108*t112*t145*t182*t353*t587*(3.0/6.25E+2)

  end subroutine WENO5_M_jab
    
end module weno_1d_mod
